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We present results from the first large-scale study of two flavor QCD using domain 



wall fermions (DWF), a chirally symmetric fermion formulation which has proven 
to be very effective in the quenched approximation. We work on lattices of size 
16 3 x 32, with a lattice cut-off of a~ x ~ 1.7 GeV, and dynamical (or sea) quark 
masses in the range m simn9e /2 

r$ ^risea r$ y^istrange' After discussing the algorithmic 
and implementation issues involved in simulating dynamical DWF, we report on 
the low-lying hadron spectrum, decay constants, static quark potential, and the 
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important kaon weak matrix element describing indirect CP violation in the Standard 
Model, Bk- In the latter case we include the effect of non-degenerate quark masses 
(m, ^ m u = m d ), finding B^(2 GeV) = 0.495 (18). 

PACS numbers: 11.15.Ha, 11.30.Rd, 12.38.Aw, 12.38.-t 12.38.Gc 
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I. INTRODUCTION 

An improved theoretical understanding of the non-perturbative aspects of QCD is in- 
creasingly important due to the continuing advance of experiments involving hadrons. Such 
understanding is an important ingredient in obtaining precise values of the parameters of 
the Standard Model and the search for new physics. It is also key to understanding the fun- 
damental properties of QCD itself which are under intense investigation at BNL, Jefferson 
Lab, FNAL, CERN, and other places. 

The fundamental theoretical tool to investigate QCD non-perturbatively is lattice QCD, 
the regularized field theory of QCD on a discrete Euclidean space-time lattice. Treatment of 
the fermion field in such a regularization has been a long-standing difficulty because flavor 
and chiral symmetry are badly broken in practical numerical simulations using conventional 
lattice fermions. A revolutionary theoretical framework to realize flavor and chiral symmetry 
on the lattice was constructed by Kaplanllf and subsequently reformulated and extended in 

q n 

two distinct ways by Narayanan and Neubergerpfl and Shamir [3J who suggested that the new 
fermions be used to study vector gauge theories, and especially QCD. These new fermions, 
known as domain wall fermions, turned out to be in a class of lattice fermions that satisfy 
the Ginsparg- Wilson relation^]. Later, Neuberger developed still another closely related 
lattice fermion called overlap fermions[5|. 

Today, both domain wall and overlap fermions are commonly used in quenched lattice 
QCD simulations, those where the fermion determinant is set to one in all path integrals used 
to calculate expectation values. In this paper, we report on the first large scale two flavor 
dynamical simulations with domain wall fermions, those where the fermion determinant is 
included in all path integrals used to calculate expectation values. From a perturbative 
point of view, this is equivalent to keeping all closed fermion loop contributions (at lowest 
order, the hadronic vacuum polarization) to all orders in all Feynman diagrams. 

To realize chiral symmetry, domain wall fermions utilize an extra fifth dimension, de- 
manding increased computational resources that are essentially linear in this extra dimen- 
sion. Nevertheless quenched QCD calculations with domain wall fermions were carried out, 
showing good chiral symmetry can be achieved with a practical number of lattice sites in 
the extra dimension (L s ~ 10) 0,0,0]. Subsequently, within the quenched approximation 
domain wall fermions were used in many calculations, showing dramatic improvement in 



scaling toward the continuum limit over conventional formulations 



0,y y 



ements 



11 



12 



14 



13JL 



were 
, and 



used in pioneering calculations of weak interaction hadronic matrix e. 
proved efficacious in non-perturbative operator renormalization lfil \v\. 

As mentioned already, large scale computations using domain wall fermions have so far 
been restricted to the quenched approximation. The problem with such calculations is 
that there is no way to systematically reduce the errors due to quenching without simply 
performing the unquenched calculation. Past experience indicates that the size of the error 
ranges between 5-10% for many observables, but can be much greater, depending on the 
observable (the critical temperature of the hadron to quark-gluon plasma phase transition 
is a well known example where the quenching error is more than 40%). In the ratio 
of pseudoscalar decay constants in the quenched approximation was found to be fx/ fn ~ 
1.13(3), different from the experimental value, ~ 1.22. On the other hand, in this work, we 
find fx/ fn = 1.175(11), as discussed in Section El Moreover, there are important physical 
objects, the flavor singlet mesons and scalar mesons, for example, which are believed to be 
very sensitive to dynamical, or sea quark effects. Ultimately, to perform accurate lattice 
QCD calculations these effects must be included. 



fermion simulations using Wilson fermions 
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Recently, large-scale dynamica 
and improved staggered fermions 21, 22] (using two light degenerate quarks and a heavier 
quark for the strange quark) have been reported. Both formulations break chiral symmetry 
and the former also breaks flavor symmetry (a severe problem, theoretically and practically). 
Nevertheless, calculations using either type of fermions, done with two or more sufficiently 
small lattice spacings, have yielded results which agree accurately with experiment for some 
observables. QCD simulations using staggered fermions are much less computationally de- 
manding, making them attractive. However, the action corresponds to four degenerate 
fermions in the continuum limit, so to simulate a single flavor requires the fourth root of the 
determinant. For this to be a legitimate procedure there must be a local action for a single 
flavor of staggered fermions which had the same determinant as this fourth root. While there 
is no proof that such a decomposition is impossible, no such action has been constructed 
to date. The flavor symmetry breaking intrinsic to staggered fermions also leads to larger 
than expected discretization errors, a problem that is significantly reduced by improving the 



naive discretization (see 
perturbation theory [3, |27 



23, 0, |3|). Even with those improvements, an extended chiral 
with many extra low energy constants corresponding to the 
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leading lattice-spacing errors must be used to do the chiral extrapolations. Without this ex- 
tra complexity, the precise agreement with experiment could not be achieved without costly 
reductions in the lattice spacing. While more computationally demanding than staggered 
fermions, the Wilson fermion formulation is still much less demanding than DWF. However, 
this action severely breaks chiral symmetry, leading to large discretisation errors (starting 
at 0(a) rather than 0(a 2 ) for un- improved Wilson fermions) and a problematic renormal- 
ization structure. In addition to this, the Wilson Dirac operator for a single flavor cannot 
be proved to have a positive determinant for positive quark mass, which is an important 
prerequisite for simulating an odd number of dynamical flavors exactly using the currently 
available algorithms. 

For these reasons, full QCD simulations with DWF are necessary as they avoid these 
theoretical uncertainties by providing a fermion formulation with both exact flavor symme- 
try, and good chiral properties. In fact , in the limit of an infinite extra dimension domain 
wall fermions posses exact chiral symmetry at non-zero lattice spacing, that is, away from 
the continuum limit. When the extent of the extra dimension (L s ) is finite, explicit chiral 
symmetry breaking is induced and is quantified in the form of a small additive quark mass 
m res . The domain wall fermion Dirac operator has the nice property that for positive mass 
and even Ls its determinant is positive, and so taking the square-root of the two flavor 
determinant is a well defined operation. This, in combination with algorithms such as the 
polynomial 28] and rational j3, 13 hybrid Monte Carlo algorithms, allows an odd number 
of dynamical flavors to be simulated exactly. In addition, if the gauge fields are sufficiently 
smooth that the underlying Wilson-Dirac operator is not in the parity broken Aoki phase 
[31 1 ) . then the theory is expected to be local. This condition has been met in quenched stud- 
ies. We refer the interested reader to 2 and references therein for a full explanation. Note 
that the same locality conditi on ap plies to overlap fermions. Thus unquenched lattice QCD 
with domain wall, or overlap jl02 |. fermions promises to be the closest non-perturbatively 
regularized theory to continuum QCD. This continuum-like (symmetry) property is central 
to the argument that despite Ginsparg- Wilson fermions being naively more expensive, supe- 
rior scaling with the lattice spacing a and greatly simplified renormalization structure may 
eventually overcome the increased computational burden since dynamical lattice simulations 
are known to scale as a - ^ 7-8 -* as the continuum, chiral, and infinite lattice volume limits are 
approached. 
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Several years ago some of us performed a study of the finite temperature QCD phase 
transition using dynamical domain wall fermions|33j]. Large explicit chiral symmetry break- 
ing was evident in, for example, the quark mass dependence of the pion mass, which is now 
believed to be caused the gauge field being sufficiently rough that there was a condensation 
of zero-modes of the four-dimensional Wilson Dirac operator appearing in the domain wall 
fermion operator (the Aoki phase); the lattice scale in these simulations was below 1 GeV. 
Here we have moved to a finer lattice spacing, a -1 ~ 1.7 GeV, and use an improved gauge 
action which significantly reduced the explicit chiral symmetry breaking in quenched calcu- 
lations by an order(s) of magnitude[12]. As discussed in Section W\ few MeV in this 
calculation, so we are confident that our simulation is not inside the Aoki phase and that, 
in fact, explicit chiral symmetry breaking is small and under precise control for a relatively 
modest L s = 12. The mass of the two dynamical domain wall quarks, m q — m/ + m res , is 
roughly in the range one half to one times the strange quark mass, meaning the residual 
quark mass due to the finite size of the fifth dimension is a small fraction of the input quark 
mass. 

This paper is organized as follows. The lattice action and its numerical implementa- 
tion with various improvements are described in Section |Hj The ensemble of gauge field 
configurations is summarized in Section IIHI Thermalization and auto-correlations in sim- 
ulation time are discussed in Section IIV1 Physical results are presented in Section IS with 
an emphasis on an analysis to next-to-leading order in chiral perturbation theory. In sec- 
tion IVII we discuss aspects of chiral symmetry related to domain wall fermion simulations 
and compare and contrast dynamical and quenched simulations. We summarize our results 
and conclusions in Section IVTTl 



II. IMPLEMENTATION OF DYNAMICAL DWF 



In this section we will discuss the algorithmic and implementation details involved in 
generating dynamical ensembles with the domain wall fermion action. As a first step we 
must precisely define the domain wall fermion action that we have used. For the domain wall 
fermion Dirac operator, I*dwf, we use the same definition and conventions as However, 
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were we simply to use the action S = S F + So with 

S F = y £ / *D DWF * (1) 

X 

and Sg representing the gauge action (in our case the DBW2 action), then we would face a 
problem: the fifth direction in the domain wall fermions framework naturally gives rise to 
unphysical heavy propagating modes in this direction, while we are only interested in the 
physics of the light mode which is bound to the domain walls. The number of such modes 
will diverge in the L s — > oo limit and so dominate the path integral. To cancel off this 
divergence we add a set of Pauli-Villars fields to the action such that S = S F + So + Spy 
with 

S PV = J2^D DWF (m f = l)<t>. (2) 

X 

Note that, while the original formulation of the domain wall fermion action included Pauli- 



34 



here we are using a slightly different form for these fields which was 



Villars fields 
introduced in 

To simulate dynamical fermions on the lattice we have a choice of many algorithms 

As in this work we are simulating an even number of dynamical 



To simulate dynamic 

00,0,0, 00. 



: lavors, it is convenient to use the well-known, exact, hybrid Monte Carlo (HMC) $ algorithm 
3J,|37J. The precise details of the HMC algorithm can be found in the cited papers, but for 
convenience we will sketch an overview here: 

• The fermionic part of the action is rewritten in terms of a set of bosonic fields using 
the relation: 

e W = det ( M ) ( 3 ) 



For this bosonic integral to converge the matrix, M, must have eigenvalues whose real 
parts are positive. Due to this condition this algorithm may not be applied to QCD- 
like theories with an odd number of flavors. However, it may be applied to theories 
where quark flavors appear in degenerate mass pairs. In this situation we may use the 
fact that 

^RD DWF R^ = D^ DWF} (4) 
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where R is the reflection operator in the fifth dimension, to rewrite det(D DWF ) 2 as 
det(D' DWF D DWF ). The matrix that will appear in the bosonic integral is therefore 

(5) 



which is the square of a Hermitian matrix and so is positive semi-definite. 

• An auxiliary field, P M , is also added to the action. This field plays the role of a 
conjugate momenta to the gauge field in a molecular dynamics evolution. 

• The evolution of the gauge field is split up into trajectories. At the beginning of 
each trajectory the pseudo-fermion field, denoted in Equation HI and the conjugate 
momenta are chosen from a heatbath. The gauge field and conjugate momenta are 
then evolved some distance in "molecular dynamics time" using a discretization of 
Hamilton's equations, which must be reversible. Several discretization techniques exist 

nhn n 

in the literature [37|, |39|, |40j ; in this work we are mainly using that introduced in |37[ , 
although we also present the results of an exploratory study of the multiple timescale 
technique of (^J. 

• If the evolution of the gauge field exactly followed Hamilton's equations then this al- 
gorithm would be complete. However, to correct for discretization errors, a Metropolis 
accept/reject step is performed at the end of every trajectory. 

Each step in the molecular dynamics evolution requires an inversion of the domain wall 
Dirac operator. While we use an even-odd preconditioned form for the operator to speed 
this up, it is, of course, still the most expensive part of the gauge field generation. In 
the following subsections we give details of our attempts to minimize both the number of 
inversions needed, and the cost of each of these inversions. 



A. New Force Term 

The initial studies of dynamical domain wall fermions [33] used separate sets of bosonic 
fields to simulate the fermionic and Pauli-Villars parts of the action. This is potentially 
wasteful, especially for a large fifth dimension, as the entire reason for including the Pauli- 
Villars term was to cancel much of the contribution from the domain wall fermion Dirac 
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operator. By using disparate sets of fields this cancellation is only apparent after the average 
over these fields; over a single trajectory the mismatch between these two pieces of the action 
may introduce large forces, and therefore large time-step errors, into the molecular dynamics 
evolution. 

In an attempt to avoid this problem we have implemented a form of the Hamiltonian for 



the molecular dynamics evolution, first suggested in |35j], that uses a single set of bosonic 
fields to estimate both the fermion and Pauli-Villars terms. To do this we simply note that 

det (DHm,)D(m,)) .,/„.■■ 1 n tmV' M 



with 



This approach needs no more memory space than the previous one and requires exactly the 
same number of Dirac operator inversions as before, although they are performed on different 
sources. However, the cancellation between the fermion and Pauli-Villars contributions to 
the molecular dynamics force now happens exactly, step-by-step in the leapfrog evolution, 
rather than stochastically. 

Table Hlgives the parameters for, and results of, a small volume head-to-head comparison 
of the old and new force terms on an 8 3 x 4 x 8 ( where the parameters are spatial volume 
times temporal length times L s ) lattice using the Wilson gauge action with j3 = 5.2, and 
Nf = 2 domain wall fermions with m/ = 0.02. As can be seen the new force term has both 
a significantly higher acceptance than the old force term when compared at the same step 
size, and a (small) reduction in the number of conjugate gradient iterations needed. The 
difference of the Hamiltonian used in the molecular dynamics evolution between the first 
and the last configuration in a trajectory, AH, is also measured. We find the theoretical 



relation to the acceptance 



39, 41 



Prob acc = erfc (^(Aff)/2) S exp (- ^S^S \ , ( 8 ) 

holds to a good accuracy. For large space-time volume, V, and small size, At, the scaling 

yf(AHY=C AH (At) 2 Vv , (9) 
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is expected, where the coefficient, C^h, is independent of volume and step size at leading 
order. By switching to the new force term, C^h is reduced by more than a factor of two, 
leading to an increased acceptance, as seen in Table HI 

A similar or somewhat better reduction of the discretization error of the new force term 
is observed for the larger lattices (16 3 x 32 x 12) and smaller couplings on which we base 
much of this work. A detailed description of the basic parameters for these ensembles is 
deferred until Section ITTT1 here we give only the details relevant for the HMC evolution which 
are summarized in Table |H] We also present the results of a short test using the old force 
term for 45 trajectories, starting from the thermalized lattice of the lightest sea quark mass. 
Using the new force term, Cah is reduced to to ~ 40% of it's value with the old force term, 
and the acceptance is increased from 56% to 77%, as shown in Table |H] An important 
observation is that C^h is almost independent of the sea quark mass for the new force term 
in our simulation. This is in contrast to the empirical assumption 
a ~ 2. 

In Section III CI we will discuss a technique that allows us to exploit this new force term 
even more succesfully. 



42 



43j, C AH oc m se °, 



B. Chronological Invertor 

For the inversion of the Dirac operator in each molecular dynamics (MD) step (the most 
ime consuming procedure in the calculation) we use chronological invertor technique of 



44j. We employ the conjugate gradient algorithm to find an approximate solution, x, of the 
inverse of the matrix D^D (here D represents the even-odd preconditioned Dirac operator), 
acting on the source vector, typ, by iteratively minimizing 

x tDWx-X^F-*U, (10) 

starting from an initial guess. In the chronological invertor technique this starting vector 
is forecast by minimizin g Eq uation HIH in the subspace spanned by the set of solutions from 
previous leapfrog steps |l03| . The precise solution is calculated by the conjugate gradient 
(CG) algorithm starting from this forecast, so the number of CG iterations is reduced. In 
this sub-section we will continue to discuss the even-odd preconditioned form of the operator. 
The first 655 trajectories in the m sea = 0.02 evolution described in Table HP were gener- 
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ated using a simple linear extrapolation of the previous two solution vectors as an initial 
guess for the conjugate gradient algorithm j^], after which we moved to the chronological 
invertor using the previous seven vectors (we found little advantage to using more than this 
number, as explained below). Figure Q shows the average and standard deviation of the 
conjugate- gradient iteration count versus the leapfrog integration step for trajectories 500 
to 655 and 3000 to 4000. The reduction in the number of CG iterations needed when using 
the chronological invertor being readily apparent. (Note: while the number of inversions 
is 52, the first and last of these are half-steps so the total distance moved in molecular 
dynamics time is 51/100). There is, however, the overhead involved in implementing the 
chronological forecast. A comparison of the time taken for producing a single trajectory 
in our particular implementation shows roughly a factor of two speed-up whereas the CG 
iteration count is reduced from that without forecasting by a factor of 2.6(1), 3.2(2), 3.3(2) 
for the m sea = 0.02, 0.03 and 0.04 evolutions, respectively. Table UTTI summarizes the number 
of CG iterations for the first 15 steps of the molecular dynamics trajectory for m sea = 0.02, 
T^sea = 0.03, and m sea = 0.04, together with the total number of CG iterations per trajectory. 
How the CG converges to the precise solution on a typical configuration of the m sea = 0.02 
enemble is illustrated in Figure 0] for various numbers of previous solutions, N p , used in the 
forecast. The precision of the forecast is improved by increasing N p . Since the improvement 
is saturated for N p > 7 for all sea quark masses used, we carried out our simulation with 
N p = 7. 

Using the chronological invertor technique, we must be careful to preserve reversibil- 
ity of the MD evolution, which is a condition for the HMC algorithm to satisfy detailed 
balance. To be precise, a trajectory starts with a gauge configuration and its conjugate 
momentum, (U^\x), Pji\x)), and evolves to another pair, (U { f\x),Pli F \x)), at the end 
of the trajectory. Starting another evolution with the latter pair with flipped momentum, 
(uji (x),— P,i (x)), the counterpart of the first configuration, (U'J; (x), — P'} (x)), is pro- 
duced. We require that U'JP (x) is the same as ujt(x) to satisfy detailed balance. 

Exact reversibility is broken in two ways in the numerical simulation. Due to round-off 
errors, the gauge links become non-unitary and are reunitarized at the end of each trajectory. 
This is not a reversible process, but only causes small changes in the link elements for the 
parameters employed in this simulation (This statement quantified below). As mentioned 
previously, a more worrying source of irreversibility is the chronological invertor. Unless 
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the convergence criteria is stringent enough so that the solution of the CG is effectively 
independent from the forecast vector, the force from the pseudo-fermion action is different 
from the corresponding flipped fermion force in the reversed trajectory, and reversibility 
breaks down. Our convergence criteria in the CG is implemented so that the relative norm 
of the preconditioned residual vector is equal to or less than a small number, R conv , 

_.. \{D\> wf Ddwf)x-$f\ < R nn 

reScG — TT— j S Rconv, (.11J 

where \ is the solution vector and $^ is the source vector. 

We define U^ N ' = U'W to be the configuration obtained by evolving with initial 
momentum P, N/2 steps followed by N/2 steps with the reversed momentum and have 
calculated the deviation between and after n steps along this path using two 
different measures: 



^• vm ) = Jrm;^ IK' w-^ML • < i2 > 

d max (U^\x)^\x)) = max | (U^\x) - Ujp(x)) | , (13) 



X,U,l,J 



where, for a generic matrix M, ||M|| 2 represents the I2 norm 45]. These are plotted in 
Figure 121 for N = 20, 40, 100, 200, 400, and 1000 using a starting configuration from the 
m sea = 0.02 ensemble (d(U^ n \ U^ 1 ') and dmaxiU^, U^) are relatively independent of the sea 
quark mass). The crucial issue is how small R CO nv must be so that the breaking of reversibility 
is negligible. In Figure El d{U^ N \U^>) and d max {U^ N \U^>) are plotted as a function of 
R C onv, for a starting configuration in the m sea = 0.02 ensemble and N = 102; the value we 
use in our simulations. For R conv > 10~ 6 the MD is not reversible. d(U^ n \ Ujp) reaches the 
edge of floating point accuracy at R con v = 10 -7 but d max (ujP\u^) is still resolved. For 
Rconv < 10~ 8 , both deviations are below single precision accuracy and comparable to the 
deviations due to reunitarization which are shown by the dotted lines in the graph. From 
these checks we choose R con v = 10 -8 as the CG convergence criterion in our simulation. 



C. Multiple Time-scale Leapfrog evolution 

In this section we discuss the multiple time-step leapfrog integration scheme of |40j . While 
this technique was not used in the main part of this work, a small study was performed to 
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test its usefulness, the results of which will be included here as they are both encouraging 
and instructive. 

In |40] a procedure was outlined for constructing leapfrog integrators for which a different 
time-step size could be used for different parts of the molecular dynamics Hamiltonian. 
The parts of this Hamiltonian which produce the dominant contribution to the leap-frog 
integration discretization error may then be treated with a finer discretization than the 
remainder. In the case where the dominant contribution to the discretization error comes 
from the gauge piece of the Hamiltonian, for which the molecular dynamics force term is 
relatively inexpensive to compute, a large improvement in the efficiency of the algorithm is 
possible. 

For the standard actions this does not seem to be the case. However, there is some 
evidence that when using the modified force term described in Section III Al the dominant 
errors are coming from the gauge part of the action. Figures 03 and |H1 contrast the individual 
contributions of the gauge, momentum, fermion and Pauli-Villars terms to the total change 
in the Hamiltonian over a trajectory for the small volume, large step size runs tabulated 
111 Table □ As can be seen, for the old force term this change is approximately the same 
size for all contributions, while for new force term the gauge and momenta contributions 
are much larger than the (combined) fermion contribution. There is also a noticeable skew 
in the distributions for the gauge and momenta pieces of the Hamiltonian, with a positive 
change in the gauge Hamiltonian being rejected much more often than a negative change. 

To test if the discretization error for the new force term is, indeed, dominated by the 
contributions from the gauge and momenta pieces of the Hamiltonian, we have performed 
a trial of the multiple gauge-step leapfrog integrator over 45 trajectories, starting from 
trajectory 1505 of the m sea = 0.04 evolution. Over this range the standard algorithm had an 
acceptance of 56%. Performing two gauge integration steps for every fermion integration step 
the acceptance moved up to 91%, confirming the gauge momenta pieces are the dominant 
cause of the discretization error. While we would like to exploit this fact by using a large 
fermion step-size combined with a small gauge step-size, we have found that in the few 
tests that we have undertaken the performance of the chronological invertor degrades as the 
fermion step-size increases by an amount that almost completely offsets the fewer number 
of inversions needed. However, this is an approach worthy of a more extended study, and 
even if it does not allow the move to larger fermion step-sizes, the gain in acceptance we 
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have observed is appreciable. 

III. SIMULATION DETAILS 

As mentionedpreviously, our simulations have been made using the standard domain wall 
Dirac operator j^J], combined with the Paulli-Villars field with the action introduced in [35]. 
For each of our three dynamical ensembles we employ two dynamical flavors, with L s = 12 
and M5 = 1.8 (we use the notation and conventions introduced in for the domain wall 
fermion action), on lattices of size 16 3 x 32. The fermion boundary conditions are periodic 
in the spatial directions and anti-periodic in the time direction. These ensembles have bare 
quark masses of 0.02, 0.03 and 0.04. The basic HMC parameters are tabulated in Table ITT1 
all Dirac matrix inversions were performed using the conjugate gradient algorithm with a 
stopping condition of 10 -8 . 

We use the DBW2 gauge action j^j with j3 = 0.80, the same action we have used in 



previous quenched studies [12J, |4j|. This action approximates the renormalization group 
trajectory by using the standard plaquette, P± x i, and a 1 x 2 rectangular plaquette, Pi X 2 : 
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-~ ^ (l-8 Cl ) TrP lxl + Cl TrP lx2 j . (14) 

This form was originally suggested by Iwasaki 0, [3, 0], who, using a perturbative 
approach, estimated a value of c\ = —0.331. In [51] a non-perturbative estimate, of 
Ci = —1.4069, was made in the quenched approximation. While we are working with 
dynamical fermions, this is the value we use; our intent being to exploit the improvement 
in the chiral properties of domain wall fermions that this gauge action provides [12j, rather 
than to be as close as possible to the renormalized trajectory. Of course, it is possible that 
the effects of the determinant will negate any such improvement. This issue is discussed 
further in section EH 



IV. THERMALIZATION AND AUTO-CORRELATIONS 



In this section we discuss the related issues of thermalization and auto-correlations for 
quantities calculated on our ensembles. Each evolution started from an ordered lattice, 
running for O(10) trajectories without the accept/reject step applied, and leaving O(600) 
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trajectories for the lattice to thermalize. The precise numbers for each evolution are given in 
Table ITVl As a check of the number of trajectories needed to thermalize the configurations 
we have calculated chiral condensate (qq)(m vau \ = m sea ) on every trajectory, using a single 
hit of a random source. Figure shows this, together with the average value for trajectory 
3000 and above, which agree well with each other by the 600th trajectory. 

For the results presented later in this work we use 94 configurations from each ensemble, 
separated by 50 HMC trajectories. For the m sea = 0.02 and m sca = 0.03 ensembles these 
configurations are taken sequentially from the thermalization point given in Table IIV1 For 
the m sea = 0.04 evolution we leave a gap of 250 trajectories, starting from trajectory 1775, 
due to a hardware error on trajectory 1772 that was not detected until after the entire 
evolution was finished. The trajectory passed the accept/reject step, and no effect was seen 
in any of the physical observables that we have calculated; nevertheless, to be cautious, we 
have allowed this gap of 250 trajectories for the evolution to re-thermalize. Table IS gives 
the results for the bare lattice values of the chiral condensate (for m va i = m sea ), and the 
r x t on-axis Wilson loop, (W(r, t)), with (r, t) = {(1, 1), (1, 2), (2, 2), (3, 3)} for these sets of 
configurations. While we do not use these values in the rest of this work, we include them 
here as they may be of utility for others trying to work at this set of parameters. 

To test this spacing of 50 trajectories we have calculated the auto-correlation function, 
defined for some observable, 0(t) (here t labels the trajectory, and runs from 1 to Nd ata ), as 

-I Ndata — i 

pit) = —— £ (0(f) - O) (0(f + t)-0) , (15) 

^ data * f—\ 

& = jj — Yl °C) » ( 16 ) 

-<» data 

on the m sea = 0.02 ensemble for the plaquette and the correlation function of the time 
component of the local axial vector current at time slice 12. The former was measured every 
trajectory, while the latter was calculated every 10 trajectories using a Coulomb gauge- 
fixed box source of size 10 and a point sink. Figures |H1 and El show both the normalized 
auto-correlation function, p(t)/p(0), and the integrated auto-correlation length, 

1 1 tmax 

T in t(t max ) = - + -^y Pi*) ' ( 17 ) 

for these two quantities, versus the separation in trajectories and the maximum separation 
over which correlations were calculated, t max , respectively. The quoted error on the in- 
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tegrated auto-correlation length is calculated using a jackknife procedure: as is standard, 
each jackknife sample is constructed by removing a contiguous group of data-points, covering 
Nuock trajectories, from the available data. On the j th jackknife sample (J 6 [0, Ndata /N b i ock )) 
we construct an estimate of the auto-correlation function, 

Pj (t) = — £ (°(0 - + 1) - O s ) (18) 



sum v j'— 1 



Ndata—t 

+ - <9J (0(f + *) - (5,) 

f=(i+l)JV Mocfc +l 

i / j N block N data , 

w „ _ w „ . I E °( f '>+ E 0(0] (w) 

f'=i t>=(j+i)N block +i 



data -<' block 



where iV s1im is the total number if terms in the two summations; this is not simply N^ata ~ 
N b i oc k - t, as if t is greater than or equal to jN block or N data - (j + l)N Uock we must drop 
the first or last summation, respectively, in Equation Estimates of the integrated auto- 
correlation length on every jackknife sample are constructed from Equation El and used to 
calculate the error in the standard fashion. Ideally, we would like to work in the regime such 
that N b i ock ^> t max , r int for a value of N b i oc k which leads to an appreciable number of jackknife 
samples. Here we use N b i ock = 100. This leads to acceptable number of jackknife samples 
(~ 50), but may be too short for a solid estimation of the error for the axial-axial data. As 
can be seen from Figures |H1 and El the integrated auto-correlation length plateaus at ~ 3 for 
the plaquette and ~ 35 for the axial- axial correlator. While these values would suggest that 
the auto-correlations may be under control for our configurations, caution should be taken 
both because their extraction is insensitive to correlations longer than O(100) trajectories, 
and because the auto-correlation length depends on the quantity. 

An important example of a quantity which displays correlations over the scale of many 
hundreds of trajectories is the topological charge. To calculate the topological charge the 
lattices are first smoothed by applying 20 steps of APE smearing with a coefficient of 0.45; 
then the topological charge is measured using a discretization of the operator 

Qtop = 32 7r 2 e 'J fc ^ r ' (20) 

which is based upon two definitions of the lattice field strength tensor, F^: that using the 
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clover leaf pattern for the simple (lxl) plaquette, 



F c 



— An 
4 



and the analogous quantity for the 2x1 rectangle, 



7 R 



An 
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where 



An (M) = - [M - M f ) 



(21) 



(22) 



(23) 



While equations ED and E21 may simply be substituted into equation EH to obtain discretized 
expressions for the topological charge (which we denote Q c and Q , respectively), they may 
also be combined such that the 0(a 2 ) errors cancel between the two definitions. This gives 
the improved topological charge operator [52! ]: 



(24) 



We ma y a lso construct a definition built up from a classically 0(a 2 ) improved field strength 
tensor |53l |: 

QT = ?0 C - 2 ~e iit ,Tr \F<F5\ + \Q R , (25) 



which is also 0(a 2 ) improved, but which has different 0(a 4 ) errors. We use this last operator 
to quote our values of the topological charge, although there would be little difference 
if we had decided to use the operator of Equation |2U for the 2565 topological charge 
measurements made on our dynamical ensembles, we observed only 3 cases in which these 
two discretizations differed when rounded to the the nearest integer. We have also compared 
our values to those calculated using the topological charge operator of [54]: for a test over 
418 measurements from the m sea = 0.03 evolution we found 90% agreement on the nearest 
integer, with the largest absolute difference being 0.76. Figure ITU1 shows these topological 
charge values versus trajectory for all of the ensembles. Note the correlations over many 
hundreds of HMC trajectories. Although the DBW2 action suppresses topology changing 
configurations (this is the reason for its improved chiral properties in conjunction with 
domain wall fermions) jl^ . considering the relatively large level of explicit chiral symmetry 
breaking observed in our simulations (see sections IVI and rVTj) . we believe the long correlation 
length observed here is due mainly to the small trajectory length HMC algorithm which is 



known to move slowly between topological sectors 
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V. PHYSICAL RESULTS 

A. hadron spectrum and decay constants 

In this work hadron masses are extracted using standard covariant fits (see, e.g. |56j |) 
of two-point correlation functions at relatively large Euclidean time interval from a single 
meson or baryon source located at t = 0, so that only a single exponential corresponding to 
the ground state is fit in each case. The source is either a Coulomb gauge-fixed wall source or 
a plain wall source. The latter corresponds to averaging over a point source on a time-slice 
after averaging over the gauge field ensemble. We use only point sinks. The wall sources 
generally have better overlap with the ground states in which we are most interested, and 
therefore lead to more accurate determination of the particle masses. However, the decay 
constants are more easily obtained from the point source matrix elements. In the case of the 
point source, we also computed the correlation functions using a source located at t = 16, i.e. 
the maximum distance from the original source on our periodic lattice, in order to improve 
our statistics. We consider zero momentum states only. All quoted statistical errors are 
estimated by fitting the correlation functions under a standard single-elimination jackknife 
procedure. 

We begin with the calculation of the residual mass, m Tes , in order to define the chiral 
limit, rrif = —m res . m res is the additive renormalization of the quark mass caused by 
mixing between the left and right handed fermions localized on opposite boundaries of the 
fifth dimension and is defined from the explicit chiral symmetry breaking term in the axial 
Ward-Takahashi identity for domain wall fermions Q]. We refer the reader to Q, Q] for 
details of the method to calculate m res . In Figure [TT] the effective residual mass, R(t), for 
n^sea — m vai — 0.02 is shown for each time slice (because the correlator is periodic, we fold 
the correlation functions about the mid-point of the lattice in the time direction which is 
why our plots range from t = to 16). It falls by about a factor of 2 over the first couple 
of times slices and then remains constant. This non-local effect of the extra dimension of 
domain wall fermions is well known 0|. Taking the error- weighted average over time 
slices 4 through 16 and linearly extrapolating the m sea = m va i points to = 0, we find 
am res = 0.001372(49) (see Figure IT2"j) . We also use the axial Ward-Takahashi identity and the 
partially conserved axial-vector current to extract the renormalization factor Za appearing 
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m We find = 0.75734(55), denned in the chiral limit, and extracted from a 

linear fit to the mass dependence of the fully dynamical points (see Figure ITS*]) . 

Next, we turn to the vector meson mass. In principle, the vector particle can decay 
hadronically into two pseudo-scalars since vacuum sea quark effects are included in this two 
flavor simulation. However, in the regime we are working the sea quarks are still relatively 
heavy, and, taking into account that two pions must have relative orbital angular momentum 
L = 1 for the decay to be allowed, it is easy to note that our vector mass is always below the 
threshold for this decay. Thus the two pion states, while present in this channel, represent 
excited states. 

Figure IT41 displays the effective mass in the vector channel, m sea = m va i, averaged over 
all three spatial polarizations. Tables IVII jVIIII give the fitted masses from the wall-point 
correlation functions for all combinations of sea and valence quark masses for a set of time 
slice ranges chosen based on the effective mass plateaus shown in Figure IT41 The masses are 
extracted from a fit to 

UmG{t) = A(e- mt + e- m(Nt - t) ) (26) 

t^oo 

= 2Ae- m(Nt/2) cosh{m{N t /2-t)). (27) 

The mass plateaus are rather poor, especially for m sea = 0.02, which is reflected in the 
values of \ 2 f° r the fits. Additional statistics are desirable, though we note our trajectory 
lengths are already quite long for dynamical fermion simulations (~ 5000). Table HXl shows 
the results of performing two fits to the data in Tables I VIII VI 1 1| a linear fit to the fully 
dynamical, m sea = m va i, points, 

M vector = a + b(m sea + m res ) , (28) 

and a partially quenched fit to all the data, 

M vector = a + b(m sea + m res ) + c{m va i + m res ) . (29) 

As can be seen from the consistancy of these two fits, the linear ansatz works well for this 
data. Considering the small mass range available, we do not attempt more sophisticated fits. 
The linear extrapolation of the three (uncorrelated) points m sea = m va i to fn = (m u + m ( j i )/2 
(later, when we discuss the pseudo-scalar meson, we will see how fa is determined) yields the 
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value of the vector meson mass corresponding to the physical p meson. From M p = 770MeV 
we find 

a" 1 = 1.691(53) GeV. (30) 

A similar analysis has been carried through for the nucleon mass. This is calcu- 
lated from Coulomb gauge-fixed wall-point two point functions of the interpolating field, 
Jn = e l i k {u lT C^$di)u k , using a positive parity projection. The results of a fit to a single 
exponential (Ae~ mt ) are shown in Table IXl Taking the negative parity projection we may 
also extract the mass of the N*, the parity partner of the nucleon, the results for which are 
tabulated in Table IXll Figure EH] shows both these quantities for the dynamical, m sea = m va i, 
points. In Figure ED we show the APE plot (Mj/M p vs. (M n /M p ) 2 where Mj is represents 
the mass of the N or N*), together with the results from quenched DWF with DBW2 gauge 
action for the nucleon. We note that the dynamical and the quenched values are in 
good agreement. Of course, in a comparison between the nucleon and rho masses, we must 
bear in mind that the rho is stable for all values of the quark masses studied here, which in- 
troduces a systematic error that could easily be 10% or more in the ratio /m p . However, 
we are encouraged to think that our error in m p may be smaller than this, since the value 
of the lattice spacing determined from the rho mass agrees quite closely with that implied 
by our calculation of f n and of the Sommer scale, r , as discussed later in this section and 



in Section IV Bl respectively. 

In Figures 1221 and we show an extrapolation to the physical point, mj — in. To 
perform this extrapolation we have used a linear ansatz for the quark mass dependence of 
the nucleon. The extrapolated values shown in the figures are taken from a fit to just the 
dynamical points, the results of which are given in Tabl dXlH As can be seen, the nucleon 
mass is two standard deviations (9%) larger than experiment. Note that our spatial lattice 
size, L ~ 1.9 fm is not small enough that we would expect to see significant finite volume 
effects in our data for the quark masses we are using; LM V ~ 4.7 for our lightest mass. 
Systematic numerical studies as well as theoretical calculations 0], on the finite 

volume effects in the nucleon mass with two flavor Wilson fermions indicate a few percent 
finite volume mass shift with similar parameters to our lightest point. However, as the p 
mass receives similar finite volume effect js?! Issj] . the ratio at the physical limit can shift as 
little as ^ 1%. Assuming that holds also for DWF, finite volume effect may be responsible for 
a minor part of the discrepancy. In addition, the systematic error associated with the linear 
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fit is at least of the scale of this discrepancy. This can be seen by comparing the diagonal 
extrapolation and two stage linear extrapolation, where the valence limit m va j — > fh is taken 
first with fixed m sea , then dynamical extrapolation m sea — > fh is performed. The result is 
shown in Table IXII| Table IXIH[ and in Figure El (triangle down) . The difference between 
the two extrapolations indicates inapplicability of the linear ansatz. This ansatz clearly does 
not properly account for observed mass dependence of the nucleon. However, the statistical 
accuracy on the nucleon mass needs to be improved before more appropriate fits can be used 
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The physical N* mass is studied similarly. The diagonal and two stage extrapolations 
yield different results, but both are consistent with experiment to within one standard 
deviation (10%) below and above. We note that the Mjy* at the lightest simulated point 
might suffer from the threshold effect, since M N + M n = 1.045(12) and M N * = 1.021(71). 
However the large statistical uncertainty prevents us from discussing further on this point. 

The analysis of the pseudo-scalar mesons will be somewhat more involved since partially- 
quenched chiral perturbation theory is at our disposal, and there are two channels, pseudo- 
scalar and the axial vector, which couple to the physical pion and kaon states. Figures HfiBTHl 
show the pseudo-scalar meson effective mass computed from the point-point and wall-point 
correlation functions. In each case we show m sea = m val ; however the cases where m sea ^ m val 
are similar. As expected, the point-point correlators exhibit plateaus that emerge at later 
times compared to the wall-point ones. As the quark mass increases, even the latter plateau 
emerges (from below) at rather late times. The statistical errors for the point-point axial- 
vector case are somewhat larger than for the others, especially at smaller valence quark mass. 
Based on the above plateaus we fit the correlation functions from i min = 9 to t max = 16. The 
results do not change by more than one standard deviation when t min is changed by two units 
in either direction, provided x 2 /dof remains acceptable which was almost always the case. 
Results for all quark mass combinations are summarized in Tables IXIVIIXVHl Fitted meson 
masses among these four methods are mostly consistent within one (statistical) standard 
deviation, except for the lightest point (m sea = 0.02, m va i = 0.01) where the mass in the 
point-point pseudo-scalar channel is almost two standard deviations higher than the rest. 
This may indicate that our statistical errors are underestimated due to the limited length 
in simulation time of our evolution (~ 5000 trajectories). 

The pseudo-scalar decay constant fp$ is obtained from the point-point correlation 
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function, either directly (axial- vector), or through the Ward-Takahashi identity (pseudo- 



scalar) |lG , 12 1 . When using a point source, A in Eq. Ellis proportional to the square of the 



decay constant. Specializing to the pseudo-scalar, 

where \ = ^7^75^ for the axial vector and ^75^ for the pseudo scalar, and Mp$ is the 
pseudo-scalar mass. We have for the lattice matrix elements 

(OlVws^) = f l P at s M PS = ^M PS , (32) 

M 2 

^ w ° fa ^^r (33) 

where the first equation defines the decay constant and the second results from it and the 
use of the Ward-Takahashi identity. The finite renormalization factor appears in the first 
equation as we use the local definition of the flavor non-singlet axial vector current which is 
not (partially-)conserved. On the other hand, no such factor appears in Eq. ESI because the 
combination (m va i + frares)(0| , ?/>75?/>|-PS') is a renormalization group invariant, protected from 
renormalization by the Ward-Takahashi identity for all a. The bare quark mass associated 
with the field ip is m va i + m res . 

Tables IX VII and fX VIII show the results for fp$, which agree well between the two methods 
except for the heaviest valence mass point at m sea = 0.02 where there is a ~ 1.5 standard 
deviation discrepancy between central values. Note that the errors from the axial-vector 
correlator are 2-3 times larger than from the pseudo-scalar, so we quote central values and 
statistical errors from the latter. 

Next, we wish to extrapolate our results to the physical light quark masses, m = 
(m u + rrid)/2 and m s (our simulation is not at a level of precision sufficient to account 
for isospin breaking effects arising from m u 7^ m^; thus we work with degenerate up and 
down valence quarks). Since we have simulated a theory where the strange quark is quenched 
and m sea 7^ m va \, the proper framework for such an extrapolation is partially-quenched chiral 
perturbation theory 66]. The next-to-leading order (NLO) partially quenched chiral pertur- 



bation theory formulae for the 



degenerate valence quarks are [62, 68]. 



pseudo-scalar mass squared and decay constant made from 



Mhd-ioop) = M 2 (l + ^) (34) 
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(35) 

6 

-[(L 5 -2L 8 )M 2 

(36) 
(37) 



(38) 
(39) 
(40) 

(41) 

The subscript U S" implies a sea quark inside the meson, A x is the chiral perturbation theory 
cut-off scale, / is the decay constant in the chiral limit, and are Gasser-Leutwyler low 
energy constants that appear in the 0(p A ) chiral lagrangian of QCD. 

We begin by fitting the pseudo-scalar meson mass to Eqs. EUandEHl To be fully consistent 
we should perform a simultaneous fit of Mps and fps as both fit functions depend on Bq 
and /. However, as we explain in the following, the fps fit is problematic. As such, for the 
the Mps fit we use a fixed value of af = 0.078, which is the result of a linear fit to the decay 
constant mass dependence. The results are summarized in Table IXIXI and, for the wall- 
point axial vector, shown in Figure ITUl We have taken A x = 1 GeV; a change in this scale is 
absorbed into the scale-dependent Lj, leaving the value of Bq unchanged. The NLO formula 
works well for quark masses up to about 0.04, where the fit quality begins to degrade. From 
Table IXVIII1 we note that a simple linear fit (lowest order in chiral perturbation theory) to 
the m = m sea = m val points, 

Mp S = c + b(m + m res ) , (42) 

is consistent with the NLO fit. This fit and data are shown in Figure [201 in fact, the 
deviations from this simple linear form are quite small. However, so is the contribution of 
the logarithm predicted by NLO chiral perturbation theory. We conclude that our data are 
consistent with NLO chiral perturbation theory for rrif ^ 0.04. 
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The NLO fit is constrained to vanish in the chiral limit, nif = — m res , as it must due to 
the universality of the low energy domain wall fermion theory Q]. Remarkably, the simple 
linear fit, which is not constrained, extrapolates to m 2 PS « (within 1 standard deviation 
of the statistical error) at m./ = — m rcs (see Fig ure l2T)|) . This is in stark contrast to the 
situation in the quenched approximation Q, lid ], where the vanishing point from a simple 
linear extrapolation occurred when rrif rs — (2 — 3) x m res . This we attributed to the presence 
of quenched chiral logarithms, ~ rn 2 ps l nm ps- Here Eq. 13*61 shows the chiral logarithms are 
much weaker, ~ m PS lnmp S . Thus this comparison of the quenched and Nf = 2 theories 
nicely confirms the predictions of chiral perturbation theory and the low-energy effective 
theory of domain wall fermions jlfll. . 

Another interesting feature of Figure is that our simulations happen to coincide with 
the region where sea quark mass effects on m 2 PS appear to be the greatest (in an absolute 
sense). Near the chiral limit, m sea = constant curves approach each other since they must 
vanish at the same place, and as m va i gets heavy, they also merge since the heavy quarks 
are insensitive to the light sea quarks. 

Now we discuss the determination of fh. fh is found from the intersection of the NLO 
fit to (amps) 2 with the line (amps) 2 = (am w ) 2 where m n is the mass of the neutral pion, 
134.9766 MeV, and the lattice spacing a is set from the vector meson mass evaluated at 
fh. This procedure is performed iteratively until convergence (in practice the number of 
iterations ^ 5). We find fh = 0.00017(11). Similarly, we find tudk = 0.0225(15) , where 
triDK is the quark mass for which a pseudo-scalar meson made of degenerate quarks has the 
same mass as the neutral kaon, m^-o = 497.672 MeV. Since the NLO formulas for M K and fx 
for non-degenerate valence quarks depend on the same parameters as in Eqs. ESI anc 
make a determination of m s (and fx) by extrapolating to this non-degenerate limit 
These equations read: 
AM 2 —1 

~Mf = N(M K - M 2 ) [ (M * " M ss^M 2 ) + (~2M K + M 2 + M 2 SS )A (M 2 33 
16 

--[(L 5 - 2L 8 )M 2 K + (L 4 - 2L 6 )NM 2 SS ] (43) 

A/V 1 , 9 n , mi + m 2 Y (m 2 qq — 2m 2 ) ( 1 , . 9 , 1 , ~ . 

-T^ = ah« 2,2 (™k - m 2 ss ) - *~ 0AT( K \ SS 2 , w) —Ao(m 2 n ) - —A (ml 3 ) 
f N16tt 2 j 2 2N(m K — m 2 r ) m 33 

N 8 
--(A (m 2 uS ) + A (m 2 sS )) + j^(L,m 2 K + L 4 Nm 2 ss ), (44) 
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Ml = 2B (m + m res ) 



(45) 
(46) 



Ml 



K — 



B (m + m s + 2m res ) 



Ml 



33 — 



2M 2 K -M 2 , 



(47) 



Note: Equat ions EH] and E31 have different logarithmic terms. The mass of the non-degenerate 
meson with one valence and both dynamical quark masses fixed to m, is also shown in 
Figure EE plotted versus the mass of the remaining valence quark. Thus, using Eq. 03] 
and the parameters determined from the degenerate formula, we find m s = 0.0446(29), our 
final value. As a consistency check, we may take this value, together with the value for m, 
and use the results of the partially quenched fit to the vector meson mass (Table UXj) to 
estimate the mass of the meson. This gives a value of M^ = 978(12)MeV; comfortably 
close to the experimental value of 1019 MeV. For comparison we have also extracted the 
strange quark mass from the linear fit. Using the three m sea = m va i points, one finds a 
value for the strange quark mass of m s = 0.04177(64), 7% smaller than the NLO value. The 
difference is easily appreciated from inspection of Figures E3 and EH and demonstrates the 
significance of the NLO analysis in this case. Finally, it is important to note that m and m s 
as quoted above are bare quark masses that correspond to the physical ir- and K- meson 
states; the renormalized quark mass is defined as Z m (m + m res ), where Z m is a scheme and 
scale dependent renormalization factor. 

We now move to the extraction of the decay constants from the point-point correlators. 
Recall that the errors on the decay constant from the axial- vector case are significantly larger 
than in the pseudo-scalar case, so we focus on the latter; our conclusions do not depend on 
this choice. In contrast to the above analysis for m 2 PS , the NLO chiral perturbation theory 
formula for the decay constant, Eqs. ESI and EH does not fit the data. The results of this fit 
are presented in Table IXXl Restricting nif < 0.03, the resulting fit is shown in Figure PHI 
While it reproduces the data used in the fit reasonably well, it misses the larger mass points 
badly. Including these points in the fit does not change the fit results significantly except 
to increase the value of x 2 . Note that the m sea = rn line bends down steeply as m va \ — > 
which yields a value for f n = 100(10) MeV, ~ 30% smaller than the physical value. 

There is just enough data to attempt a restricted NNLO fit, in which we include all of 
the 0(p 6 ) analytic terms: Cim 2 ea , C 2 m 2 al , and Ci 2 m sea m va i. This fit is not a systematic ap- 
plication of chiral perturbation theory, as we do not include the (un-calculated) logarithmic 
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terms which also appear at this order. However, performing this fit allows us to investigate 
the utility of moving to the next order, and estimate the size of the terms needed. The 
results of the NNLO fit are summarized in Table IXXll While the value of y 2 is acceptable, 
the errors on the fitted parameters are extremely large, especially in the case of G\. 

The basic problem is that the coefficient of the log term, which has been calculated 
analytically in the continuum (Eq. IH7|) . is inconsistent with our data, which appears to be 
mostly linear. The additional NNLO terms in the fit do act to reduce the effect of this log, 
but this approach is both badly constrained and non- systematic. One interpretation of our 
data is that the quark masses we are using are sufficiently heavy that the NNLO terms are 
important relative to the NLO terms. In this case, the only way to do better is to simulate at 
even lighter values of the quark mass where (presumably) chiral perturbation theory works 
well. 

There are several other interpretations: both lattice artifacts (0(a 2 ) in this study) and 
finite volume effects modify the coefficient of the chiral log, possibly making it smaller. In 
the former case we expect this effect to be of the order of a few percent, while including finite 
size effects such as the ones in 69] should not change the fits significantly since our smallest 
mass still corresponds to m ppL ^ 3. Neither of these effects therefore seem large enough to 
explain the discrepancy. In 7fJ, it is also suggested that using a physical parameter, such 
as fx, as the chiral coupling rather than /, leads to a better behaved chiral expansion. This 
approach will have the effect of significantly reducing the coefficient of the chiral logarithm 
when applied to our data. If we allow the coefficient of the continuum chiral log to be free 
parameter in the NLO equation, then we are able to make a good fit (x 2 /dof = 0.48(38)). 
However, we find that this coefficient is 0.2(4), instead of the value 1 predicted by the 
continuum theory, which is a large deviation. 

A remaining alternative is to forsake most of the higher order terms and do a simple 
linear fit. Again, this is not systematic, as this leaves out the (known) NLO logarithmic 
term. We fit the data with three independent terms, /, and terms proportional to m sea and 
m va i (in other words, we simply set the coefficient of the log term in the NLO formula to 
zero): 

m 1 + m 2 + 2m res 

Ips = a J + ci h c 2 (m sea + m res ) (48) 

As mentioned before, such a fit - the results of which we summarize in Table IXXIII and 
Figure ED - actually describes the data quite well. For comparison, / is also extracted 
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from a simple two parameter linear fit to the m sea = m va i data points (see Table IXXIII and 
Figure l2*2*|); these two linear fits agree well. Using the three parameter linear fit, restricted 
to m val < 0.04, we find f n = 134.0(42) MeV, f K = 157.4(38) MeV, and f K /f w = 1.175(11), 
while the Particle Data Book gives / w = 130.7 MeV, f K = 160 MeV, and f K /U = 1-224. 

While the inability to fit our data to the predicted NLO chiral perturbation theory form is 
discouraging, it is not unprecedented: this problem also exists for the current Wilson fermion 
[rsl \v\ . and staggered fermion [3] simulations. The former simulations, for dynamical 
masses in the range such that M^/Mp = 0.6 — 0.8, saw little evidence for the existence of the 
chiral log, extracting final numbers using a polynomial ansatz for the mass dependence ( [zj] 
discusses the matter further). The latter simulations, which include much lighter dynamical 
masses, also cannot fit the NLO formula for comparable quark masses. They quote final 
results from a restricted NNLO fit which, as above, leaves out the logarithmic terms that 
should appear at that order. The advantage of this work, however, is that we have a much 
cleaner extraction of the decay constants due to both flavour and chiral symmetry being 
respected at finite lattice spacing, allowing us to be confident that this discrepancy from 
continuum NLO chiral perturbation theory is a physical effect. 

B. Static quark potential 

In this section we discuss the extraction of the static quark potential, V(r), from the 
Wilson loop, W(f, t). To improve the signal we smear the operator in the spatial coordinates 
using n smear applications of APE smearing [73]. 

U[{x) <- Proj SU (3) [Ui(x) + c smear U\ staple \x)\ (i = 1, 2, 3) , (49) 

where (x) is the sum of four spatial staples. Using these smeared links we con- 

struct the spatial path between the infinitely heavy quark and anti-quark by employing the 
Bresenham algorithm 

To cut down on short-term noise, we measure the static quark potential more frequently 
than every 50 trajectories, although our final results are obtained from a block-average 
over 50 trajectories. Fig. [2H shows the integrated autocorrelation times of a selection of 
Wilson loops, which are < 25 trajectories for r < 8. In total we use 941, 559 and 473 
configurations for m sea = 0.02, 0.03 and 0.04 respectively; these configurations are separated 
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by five trajectories for m sea = 0.02, while one in every ten trajectories is measured for 
m sea = 0.03 and 0.04. Trajectories from 1775th to 2025th of m sca = 0.04 are abandoned due 
to the hardware error described in Section IIV1 
We use the theoretical formula 

(W(r,t)) 



C(f) exp[-V(r)t] 



(W(r,0)) 

which should hold for large t, and (following |l2j|) extract V(f, t) from 

(W(r,t)) 



together with C(f) from 



V(r) = lo{= 



(W(r,t+1)) 
(W(r,t)) t+1 



(50) 



(51) 



(52) 



(W(r,0))(W(r,t+l)y ' 
To maximize the overlap factor, C(r), we have explored the two dimensional parameter 



space for the smearing, (c 



smear j ^ smear ) > 



in the range < c s 



< 1, < n s 



< 60. We 



conclude (c smear ,n smear ) = (0.5,20) is a reasonable choice for all of r we use. 



The dependence of the potential to the temporal length of the Wilson loop is carefully ex- 
amined to control the contamination from excited states. We found V(r) and C(r) extracted 
from t > 5 and r < 8 is relatively independent of t, and we therefore extract our results from 
this range. This can be clearly seen in Figure EEl which shows V(r) for m sea = 0.02. The 



[ which shows V(r) for m se 
effect of the positivity violation in improved gauge actions 7^] for small t is evident in the 
graph: V(r) extracted from t < 3 approaches its asymptotic value from below for r = \/2. 
This was also observed in quenched simulations 17151 . and the effect on C(r) is discussed in 

n n 

|76L 1771 ] : our selection of temporal length is set as large as possible to exclude this lattice 
artifact. The potential extracted from t — 4, 5 and 6 is shown in Fig. EBlfor the m sea = 0.02 
ensemble versus r. We see no evidence of string-breaking for large r. 

For ease of comparison with other work, we also fit the potential to the form 



V{r)=V cont {r)-lSV(r), r=\r\ 



cont \ 



5V(r) 



a 

V + or 

r 

I 

r 



where [l/r\ is the lattice Coulomb potential 



78, 



79 



80] 



" dk 3 



exp(ik ■ r) 



n 8tt2 £. sin 2 (^/2) - 4 Cl sm 4 (fc 4 /2) 



(53) 
(54) 

(55) 



(56) 
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This form is often used in an attempt to correct the breaking of the rotational symmetry 
of the lattice. The correction term, SV(r), describes the corresponding data, |V( data )(r) — 
V con t(r)]/l, qualitatively well as seen in Figure U7\ in which V cont {r) and I are from the fit 
result using t = 5, r G [a/3, 8] on m sea = 0.02 ensemble. We also note that the fit parameters, 
especially a, become less sensitive to the selection of r min by adding the lattice Coulomb 
term, although the errors become larger. 

it results (both with and without the lattice Coulomb term) , together with Sommer 



The 
scale 



/ 1.65 — a 

n> = \/— ^— , (57) 



are presented in Table. IXXIIII As can be seen in Figure EH1 a major source of systematic 
error is the selection of t, with the r max dependence being almost negligible. To be precise, 
one can see a mild but continuous increase of tq as t becomes large. At t — 7, the signal 
to noise ratio becomes poor, and the statistical error of the data at t — 6 is less controlled 
as seen in Figure We therefore choose to take our central values from t = 5 and 
r £ [ r 'min, r max] = [\/3, 8], quoting a systematic error due to the selection of t, as well as fit 
range by the shift of central values of these parameters. The selection is varied in either 
direction in t,r min ,r max at once. The fit ranges r min G [a/2, a/6], r max G [7,9], and t = 5, 6 
are swept. More detailed results including the direct evaluation of the force, VV(r), and 

nn 

comparison with quenched simulation will be presented elsewhere |76J,|82J. 

Although in this paper we use the hadronic observables to set the lattice scale, we may 
also set the scale from tq. The sea quark mass dependency of r$ is shown in Figure OH 
This mass dependence is consistent with that observed using other fermion formulations 

q y a Q a n,^ e „ ing * ^ ^ we obtaui . ^ * ro m 

lattice unit of 

r | m __ m _ = 4.278(54) . (58) 

This value is obtained from the fit without the lattice Coulomb term and the systematic 
error in the second parenthesis is estimated by the various choice of fit ranges. We note that 
the large positive shift of the central value (+174) is due to the rise of the r at r min = a/3 
from t = 5 to t — 6. This could be a sign of the remaining excited states contamination but 
it is less conclusive with the comparably large statistical error at t = 6 as seen in Figure I2TJ1 
We note that vq depends on m sea so mildly that linear fit of 1/ro yields a very close value, 
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r o|m sca ^-m rcs — 4.287(58), at the chiral limit. The fit with the lattice Coulomb term yields 
similar value, ro = 4.235(56). All of these three numbers are within the quoted statistical 
error. 

For our final results we use the continuum fitting form, and a linear fit to the mass 
dependence; taking r = 0.5 fm we get 

a" 1 = 1.688(21) GeV , (59) 

n 

which, in contrast to the situation in the quenched approximation |12|, is consistent with the 
value, a" 1 = 1.69 1(53) GeV, extracted from the rho meson mass. 



C. Kaon B parameter 

When studying the mixing of neutral kaon in the standard model, it is necessary to 
calculate the low energy matrix element, in QCD, of the effective weak interaction operator 



L l = sj^l - 7 5 )ds7 M (l - 7 5 )d, (60) 
between kaon states. This is usually quantified in terms of the kaon B-parameter, 

B = (K°\si»(l-l 5 )dsMl-l 5 )d\K°) 
K ~ l(K°\A 4 \0)(0\A,\K°) ' 1 ) 

(^1^(1-75)^^(1-75)^1^°) 

3 JK m K 

Before describing the lattice calculation of Bx, we address the (unphysical) mixing of 
Oll with wrong chirality dimension six operators. Such mixings are allowed in the chiral 
limit due to the explicit breaking of chiral symmetry of domain wall fermions with finite 
L s . However, since such explicit breaking is small, it is useful to understand the order of 
magnitude of these mixings in terms of am res . This allows us to argue they may be neglected 
in our calculation. The framework for understanding this mixing within a low energy effective 
theory of domain wall fermions was laid out in , but there the explicit example of B% 

was not discussed. We begin by adding a spurion field Q to the low energy effective action for 
domain wall fermion QCD, whose only effect is to modify the symmetry breaking terms in 
the action so they transform in the same way as a conventional mass term. Each spurion field 
in the modified action or effective weak operator carries a factor of am res . After determining 
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the dependence of the correlation function on f2, and therefore am res , we set — > 1 to 
recover the low energy theory of domain wall fermions. In the exact theory, the original 
four-quark operator transforms as a (27,1) dimensional representation of SUi(3) x SUr(3) 
chiral symmetry which must also be true for the modified wrong chirality operators if they 
are to mix. All such operators [86] have at least two right-handed quark fields, so at least 
two spurion fields are needed to formally rotate these into left-handed fields. Thus, to lowest 
order the mixing is ~ (am res ) 2 , or (9(1CT 6 ); this is small enough that it may be neglected in 
the present calculation. We note that the leading explicit chiral symmetry breaking in the 
matrix element is still (9(am res ), however. This is just the statement that the chiral limit for 
all low energy observables is m/ = — m res [10]- Once this trivial shift is taken into account, 
the error on Bk is again 0((am res ) 2 ). 

Details of the method used to calculate B K can be found in [l^. We employ the so- 
called conventional method (Eq. 161(1. since difficulties of the quenched approximation, such 
as significant contamination by topological near- zero modes, do not apply Q|. The lattice 
B parameter for arbitrary quark mass, denoted Bps-, is tabulated in Tables IXXI VI and IXX VI 
IXXVII( for degenerate and non-degenerate valence quark masses respectively. The values 
shown are averages over the time slice range for the operator insertion, t op = 14—17. The 
source and sink pseudo-scalar meson time slices were fixed to t = 4 and 28, respectively. 
These values were chosen based on the quenched calculation in [15] (a -1 f» 2 GeV) where 
they lead to reasonable plateaus, which is also the case here (we also averaged the value of 
Bps over a larger time slice range, t op = 10 — 22, and found good agreement in all cases). 
This can be seen in Figures l3~Tll3*3l Note however, that there are large fluctuations in the 
plateaus for the smallest valence quark masses. 

To extract B^, we fit our data for Bps to the predictions of NLO chiral perturbation 
theory, and extrapolate/interpolate to the physical quark masses. Our preferred way to do 
this is, of course, to take the (non-degenerate ) limit 

m sea -> m ; m^i m s ; m va i i2 -> m, (63) 

corresponding to dynamical, degenerate up and down quarks, and a quenched strange quark. 
However, as previous work has only used degenerate valence quark masses, extrapolating to 
the point 

mvai,! = ra va i,2 = m DK , (64) 
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we also present results from this method, so we may compare the two techniques. 

The degenerate valence quark mass data is plotted in Figure OH As can be seen, Bps 
displays a rather weak, but noticeable, dependence on the sea quark mass. To be precise: 
as m sea decreases so does Bps- The predicted dependence of Bps on m sea and m va i, when 

Bps = bo(l- —^r U M 2 log ¥P\)+ (6i - h)M 2 + b 2 M 2 ss . (65) 



(4vr/) 2 V " A 2 

Fits using Eq. IHElare summarized in Table IXX VI I II Eq. ESdoes not simultaneously fit all of 
the data very well. This is evident from FigureEU the results of the fit for 0.02 < m va i < 0.04 
are shown, and under-predicts both the lightest and heaviest points. However, as noted 
previously, we see large fluctuations in the plateau for the lightest valence quark masses, 
and the systematic error on these points probably outweighs the quoted statistical error. 
Therefore, for our final results we exclude both the m va i = 0.01 and m va i = 0.05 points; the 
latter in an attempt to stay in the region for which NLO chiral perturbation theory is valid. 

We now move to the extraction of Bk including the effect of non-degenerate valence 
quarks. In this case sea quark dependent log terms appear, as well as many new non- 
degenerate valence quark terms: 



M? 9 ( , „ , M? 9 , ,,, . 9x , 1 + e 
_^_ 2(3 + £2)log _. 



Bps = Ml + 77TlM-2(3 + e 2 )log-f - (2 + e 2 ) log (1 - e 2 ) - 3e log ) (66) 



2 1 I 3 /2-e 2 1 + e 
^WJf\N MsS \— l0t —e 

4 M?2 U + t-*-* log l±f + ^ lor MUl ~ e] 



N V e 1 - e A 2 

M 2 

+h M 2 2 + b 2 M 2 SS + h M 2 PS (-2 + -^f) . 

In the above 

e ss m2 ~ mi , (67) 
mi + m 2 + 2m res 

M 2 PS = 2B (m l + m Ies ), (68) 
M 2 2 = 2B — \ -, (69) 

where mi, m 2 (mi < m 2 ) are the valence quark masses. The log terms were computed in 
js7|. Bps is shown in Figure EH for both the degenerate and non-degenerate mass points, 
plotted as a function of (mi + m 2 )/2. As can be seen, the non-degenerate points appear to lie 
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on a smooth line connecting the degenerate mass points, a clear sign that the non-degenerate 
mass effects are small. Fits to Eq. EH1 are summarized in Table IXXIXI Note that the value 
of Bk is determined from the fit to Eq. IHE1 evaluated at m s and rh. As in the degenerate 
case, we use valence quark masses in the range 0.02 < m va i < 0.04 for the extraction of our 
final result. 



It is customary to quote Br in the NDR-MS 1 scheme at /i = 2 GeV. We accomplish this 
in two steps. First we renormalize non-perturbatively in the RI scheme at a scale \i ~ 1/a 
(Note that this extraction makes use of the pertubative two-loop, continuum, results 



for the running ), and then match to MS in the continuum using one-loop perturbation 



theory 



891 ]. Details of our method are given in 



47, tea. With Z 



MS 



(adding errors in quadrature) 



for the degenerate case, and 



5^(2 GeV) = 0.509(18). 



B A K IS {2GeV) = 0.495(18), 



0.93(2) 88J we find 



(70) 



(71) 



for the non-degenerate case. We take this as our final value for Bk- While these two 
numbers agree within the quoted statistical error, taking the jackknife difference shows that 
the 3% difference between the non-degenerate and degenerate extraction is statistically well 
resolved. 

y 10% smaller than recent 
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15 



47 



901 ] and roughly 



This two flavor, non-degenerate valence quark result is rough 
quenched values obtained with domain wall or overlap fermions 

20% Sm a„er than the inched va, U es reported in Q Q <Ko gU t-S US8ki »d and over.ap 
fermions, respectively). We caution that there are as yet unquantified systematic errors in 
our determination of Bk : most notably non-zero lattice spacing and finite volume effects. 



Recent quenched studies 



11 



471 ] find these to be on the order of 5% each, though they could 



differ in the dynamical case. 

Of course, the quantity of direct relevance to experiment is B^ at the kaon mass, as 
quoted. However, in passing, we also mention our value for this parameter in the chiral 



limit, B K (MS, 2GeV), as it provides a useful comparison with phenomenological models 

where calculations in the chiral limit are often under better control. We find, 
B K (MS, 2GeV) = 0.241(10). This should be compared with a value of 0.267(14), obtained 
in our previous quenched study jis| . 
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In summary: our study shows that the inclusion of sea quarks and non-degenerate valence 
quarks tends to lower the value of B K by roughly 10% and 3%, respectively, which represents 
a very important step in estimating systematic uncertainties in Bk- We note that some time 
ago Nt = 2 staggered fermion calculations found no sea quark effects on Bk outside of quoted 
errors |9fil I97I . These studies used unimproved staggered fermions which have large lattice 
spacing errors; thus we consider our new determination to be more reliable. It should also 
be noted that early Wilson fermion resuto Q, as weU as the more recent l99f suggest a 
small decrease in Bk when including the effects of dynamical quarks. In |l00 | the effect 
of sea quarks on Bk was estimated in chiral perturbation theory from the difference of the 
chiral logs between the quenched and dynamical theories. This comparison suggested that 
the sea quark effects increase the value of Bk, however it assumed that the analytic terms 
remain the same between the quenched and dynamical theories and so, again, we consider 
our determination to be more reliable. 



VI. CHIRAL SYMMETRY 

In this section we will discuss our understanding of the (relatively large) breaking of chiral 
symmetry observed in this work, starting with our choice of gauge action. As mentioned pre- 
viously, the DBW2 gauge action has been used successfully in the quenched approximation 
to improve the chiral properties of domain wall fermions. To give an example: when com- 
paring the Wilson, Iwasaki, and DBW2 actions for inverse lattice spacings of approximately 
2GeV, the residual masses for L s = 16 are around 3MeV, 0.3MeV and 0.03MeV respectively 



12( . This dramatic improvement can be ascribed to the interplay of the positively weighted 
plaquette term and the negatively weighted rectangle term leading to a strong suppression of 
dislocations of the lattice. While we would like to exploit this mechanism in the dynamical 
case, it is not a priori obvious what form of bare gauge action we should take such that, when 
combined with the effects of the fermion determinant, we have an effective short- distance 
gauge action with a similar form to that used in the quenched approximation. Over large 
(physical) distance scales it is well-known that the effect of adding the determinant is to 
smooth out the gauge field, and so, in order to perform a dynamical simulation with the 
same lattice spacing as a quenched simulation, we must increase the gauge coupling. The 
question we are interested in here, however, is how the inclusion of the fermion determinant 
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modifies the short- distance properties of the gauge field. One particular worry, for example, 
is that the addition of the fermion determinant may effectively induce a rectangle term of the 
opposite sign to the one in the gauge action, leading to reduced suppression of dislocations. 

To study this we solve for the short distance effective gauge action using the Schwinger- 
Dyson equations, following the approach of j^jj. Briefly summarizing: for an ansatz of the 
effective gauge action of 

S = /3 a S a , (72) 

a set of operators, Oi, and some variation of the gauge fields, 5, the Schwinger-Dyson 
equations read 

(<50 i ) = -(O i [^ a ])/3 Q - (73) 

Calculating {50 j) and (Oj [<55' a ]) for the same number of independent operators as there are 
independent terms in the ansatz for the gauge action, this equation may be solved to give 
(3 a . In the following we use 

6 i 

where U\ denotes a specific link and Gf the sum of the (forward) staples of type a for link 
I, and the variation 

SfUi = -ie a X a U t (75) 

for observables 

Of = ImTr [\ a UiG]\ , (76) 

and then sum Eq. 1731 over a and I. 

Some idea of the utility of this method can be gained by studying quenched configurations. 
Obviously, should the ansatz of the gauge action include the terms that constitute the 
action used to generate the ensemble, then this is precisely the result this method will 
give. However, in Table IXXXI we show the results of applying this method using a quenched 
ensemble of 404 configurations of size 16 3 x 32 generated with the DBW2 action with (3 = 1.04 
- corresponding to an inverse lattice spacing of ~ 2GeV - and solving for the effective 
plaquette gauge action. Using an observable based on the simple plaquette staple we get 
an answer of (3p ~ 9. An ensemble generated with this plaquette gauge action would have 
a much finer lattice spacing than the one we are studying; our interpretation of this result 
is that, at the short distances probed by this observable, lattices generated with the DBW2 
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action are much smoother than those of the same lattice spacing generated with the plaquette 
action. As can also be seen, as the size of the observable used to calculate the effective f3p 
is increased to the 2x1 planar rectangle; 2x2 square and then 3x3 square, the calculated 
value of (3p decreases towards, presumably, the value that would generate one ensemble with 
the same lattice spacing as the one we are studying (/3 « 6). 

Table lXXXII shows the results of applying this method to each of our dynamical ensembles, 
using relatively localized observables. In each case we see that the short distance properties 
of the gauge field are dominated by the bare form of the gauge action: the effective gauge 
action, up to some very small deviations, is identical to this bare gauge action. We may 
therefore qualitatively understand the larger observed chiral symmetry breaking for the 
dynamical lattices versus quenched lattices as being due to the smaller value of the (bare) 
(3 used in the gauge action leading to reduced suppression of dislocations; the effects of the 
determinant seem to only become significant for larger distance observables. 

While the chiral symmetry breaking observed for the ensemble studied here is small 
enough to be a negligible correction to the physical quantities we are calculating, it is still 
interesting to know how the residual mass varies both with the size of the fifth dimension 
and the gauge coupling. A reasonable estimate of the former may be made by calculating 
the residual mass for different valence values of L s . Table IXXXIIll shows this for fifth 
dimension of extent 8 to 32 for the m sea = 0.04 ensemble at the dynamical point. For 
L s ^ 12 the residual mass is calculated on 45 configurations, leaving 100 trajectories between 
configurations starting from trajectory 1005 (we skip trajectories 1805 and 1905 due to the 
hardware error mentioned in Section HVjl . Figure 137)1 shows these dynamical results, and 
compares them with the residual mass on three quenched ensembles: Wilson (3 = 6.0 [10], 



Iwasaki (3 = 2.6 M and DBW2 (3 = 1.04 



12|; all of which have inverse lattice spacings 



~ 2GeV. As can be seen, while the residual mass in the dynamical simulation is much 
larger than that in the comparable quenched DBW2 calculation, it is still smaller than the 
value when using the Wilson gauge action in the quenched approximation. 

We have also made a few explorato ry st udies using stronger couplings for the DBW2 



gauge action in dynamical simulations 



lOlJ, namely (3 = 0.75 and 0.70. For both these 



simulations we used M 5 = 1.8, L s = 12 on 16 3 x 32 lattices, as with the rest of this work. 
Table IXXXIII gives values of the input quark mass, number of configurations collected, the 
rho meson mass, and the value of the residual mass. As for each value of the coupling we 
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used only one value for the dynamical quark mass, the latter two quantities are quoted in 
the valence chiral limit. While the residual mass at (3 — 0.70 is not prohibitively large, 
caution should be taken; we may also look at at the spectral flow of the Hermitian Wilson 
Dirac operator. A transfer matrix in the fifth dimension for domain wall fermions may be 
written in terms of this operator, with the success of the domain wall fermion mechanism 
being dependent on the existence of a gap in the spectral flow for negative Wilson masses. 
Figure shows a typical spectral flow for the (3 = 0.8 ensembles, while Figure EH shows 
typical spectral flows for all the gauge couplings we have studied. As can be seen, the gap 
in the spectral flow rapidly closes as we move to stronger gauge coupling. This is, perhaps, 
the main challenge for the domain wall fermion/overlap approach: the computational cost 
of lattice generation falls so quickly as the lattice spacing increases that the extra cost of 
domain wall fermions versus other approaches could easily be amortized by working on 
slightly coarser lattices. However, for the current formulation, the domain wall mechanism 
begins to fail as we move to lattice spacings significantly coarser then a~ l ~ 2GeV. There 
have been some preliminary attempts to modify the gauge and/or fe rmion action in such 
a way that the domain wall mechanism persists at stronger couplings 10ll |. These have so 
far met with limited success, but this is clearly a direction of research which has not been 
exhausted. 



VII. CONCLUSIONS 

We have presented a large scale lattice QCD calculation using two flavors of dynamical 
domain wall fermions with small quark masses on lattices with large volumes. Domain wall 
fermions possess exact chiral symmetry in the limit L s — > oo even at finite lattice spacing 
- a symmetry fundamental to much of the physics of QCD. The ~ 3 x 5000 Monte Carlo 
trajectories in this work obtained with light dynamical quarks and small residual quark mass, 
m strange /2 ^ m sea + m rcs ^ ^strange, represent a substantial computational undertaking that 
took almost two years to complete. 

The first physics results to come from this endeavor are both interesting and encouraging. 
The decay constants f n = 134.0(42) and fx = 157.4(38) agree with experiment well within ^ 
5% statistical errors. Their ratio, determined to 1% statistical accuracy, fx/ fn = 1.175(11), 
slightly under-predicts the experimental value [1.223] obtained from the Particle Data Book. 
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The Kaon B parameter is the hadronic matrix element of the AS = 2 weak interaction 
operator that governs neutral kaon mixing. It is required to determine the level of indirect 
CP violation in the kaon system that is predicted in the Standard Model. We find that 
inclusion of sea quarks tends to decrease the value of Bk relative to our comparable quenched 
calculations H] by about two (statistical) standard deviations, or roughly 10%. Non- 
degenerate valence quark effects further lower the value by 3%. At a renormalization scale 
of 2 GeV in the continuum MS scheme, we find B^ s = 0.495(18) (statistical error only) 
which is significantly lower than previously reported quenched values and could impact the 
extraction of the phase of Standard Model CKM quark mixing matrix if this value accurately 
describes the chiral, infinite volume and continuum limits. 

Besides physical results, two important and closely related theories, the low-energy ef- 
fective theory of domain wall fermions and (partially-) quenched chiral perturbation theory, 
are strongly supported by the pattern of explicit chiral symmetry breaking in our calcula- 
tions. In quenched calculations the pion mass squared, when fit to a simple linear function 
of the quark mass, significantly over-shot the chiral limit, rrif — —m res . Through careful 
study, it was shown that this effect could be explained by a prediction of quenched chiral 
perturbation theory, that is to a logarithmic singularity unique to the chiral approximation. 
When accounting for this term, was shown to vanish at the correct chiral limit. In the 
two-flavor case this offending logarithm does not appear. Such a logarithm appears only at 
higher order (see Eq. l3*Uj) and does not effect the chiral limit. The net result is that a simple 
linear extrapolation of should come closer to the true chiral limit. That this indeed 
happens was demonstrated in Figure mil and Table IXlXl 

These initial results, obtained from two flavor calculations with three relatively heavy 
sea quark masses, on a single volume and lattice spacing, and fifth dimension L s = 12, 
while encouraging, may still suffer significant systematic uncertainties. We stress that the 
use of next-to-leading order partially-quenched chiral perturbation theory was crucial in 
our analysis of and Bk, where it worked reasonably well. The analysis for the decay 
constant was more problematic. This lead us to quote results for f v and fx from linear fits, 
i.e. the NLO analytic terms were included but not the logarithms. In addition, while the 
number of trajectories in our study is large from a historical perspective, large scale auto- 
correlations which can only be detected in longer runs may still be present. In fact, here 
as in all dynamical simulations, the number of trajectories studied is determined at least 
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as much by the amount of available computer resources as by established scientic criteria. 
All of these issues must be further addressed by future calculations. The proven scaling 
behavior of domain wall fermions gives us confidence that the results presented here provide 
a solid foundation on which to build. 
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Tables 

TABLE I: Small lattice comparison of HMC evolutions. All these evolutions use the Wilson gauge 
action with (3 = 5.2 and two flavours of Domain Wall Fermions with a bare mass of m sea = 0.02. 
Force Term At Steps/Trajectory Trajectories Acceptance CG- iterations/Trajectory Cah 
Old 1/64 33 1000-1880 87% 8336 26.5 

Old 1/32 17 1000-1929 59% 4310 30.0 

New 1/32 17 1000-1936 79% 4179 12.9 



TABLE II: Parameters for the large lattice HMC evolutions. The averaged elapse time per one 
trajectory of our paricular implementation on 32 mother boards (32MB ~ 100GFLOPS theo- 
retical peak speed) or 64 mother boards (64MB ~ 200GFLOPS) QCDSP, and the observed 
acceptance in the Metropolis test are also quoted. The CPU time includes the computation 
time for the chiral condensation, (qq)^ =m , and the r x t on-axis Wilson loop, (W(r, t)} with 
(r, t) = {(1, 1), (1, 2), (2, 1)}, but does not include the time for I/O. The scaled squared energy dif- 



ference between the first and the last configuration in a trajectory, Cah = \J ({AH) 2 ) /V /(At) 2 , 
is quoted with the standard deviation error. "0.02 (OLD)" are the results of a small number of 
trajectories using the old force term. 



m sea At 


Steps/Trajectory Trajectories Acceptance 


Cah 


Time/Trajectory (machine) 


0.02 (OLD) 1/100 


51 


45 


56% 


39(4) 




0.02 1/100 


51 


656-5361 


77% 


16.2(2) 


0.8784(6) hours (64MB) 


0.03 1/100 


51 


615-6195 


78% 


15.8(1) 


0.8324(4) hours (32MB) 


0.04 1/80 


41 


625-5605 


68% 


16.4(2) 


0.7116(2) hours (32MB) 



45 



TABLE III: Average cg-count and standard deviation (shown in square brackets) versus molecular 
dynamics step for trajectories 3000 to 4000 of each ensemble, together with the average for the 
total number of cg-iterations per trajectory. 



Leapfrog 


step m sea = 0.02 


m sea = 0.03 


m sca = 0.04 


1/2 step 


715[11] 


513.9[49] 


401.5[32] 


step 1 


626[11] 


435.0[49] 


340.4[32] 


step 2 


543[11] 


363.3[47] 


284.7[29] 


step 3 


477[13] 


297.6[50] 


231.7[32] 


step 4 


411 [19] 


227.3[90] 


172. 3 [48] 


step 5 


346[15] 


181.1[61] 


137.2 [46] 


step 6 


280[12] 


175.0 [54] 


128.2[41] 


step 7 


277[13] 


157.9[58] 


120.5[29] 


step 8 


269[10] 


149.7[39] 


117.8[21] 


step 9 


274[14] 


146.0[60] 


116.4[48] 


step 10 


275[15] 


154.5 [74] 


115.6[51] 


step 11 


275[14] 


158.3[47] 


127.3[28] 


step 12 


269[17] 


153.6[78] 


122.0[65] 


step 13 


282[15] 


152. 6 [76] 


120.5[52] 


step 14 


279 [17] 


158.8[76] 


122.1[51] 


Total 


16014[396] 


9214[130] 


5964[71] 



TABLE IV: Evolution details and number of configuration used for physical results. 
m sea Until Accept /Reject Thermalisation Trajectories Configurations 
0.02 28 656 5361 94 

0.03 25 615 6195 94 

0.04 5 625 5605 94 
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TABLE V: The bare chiral condensate, (QQ)m va , 1 =m se& ^ and the rxi on-axis Wilson loop, (W(r,t)}. 
m sea N conftgs (qq)j val=maea (W(l,l)} (W(2,l)} (W(2,2)) (W(3,3)) 
0.02 94 0.002542(11) 0.646706(50) 0.407777(80) 0.17601(10) 0.033633(80) 

0.03 94 0.0034419(98) 0.646594(47) 0.407629(82) 0.17572(10) 0.033287(81) 

0.04 94 0.0043255(81) 0.646561(48) 0.407562(71) 0.175615(97) 0.033330(81) 



TABLE VI: The fitted vector meson mass from the wall-point correlation function for m sea = 0.02. 



m sea 


m va i 


fit ran£ 


;e x 2 / dof 


mass 


0.02 


0.01 


5- 


16 


1.9(10) 


0.511(10) 


0.02 


0.015 


5- 


16 


2.0(11) 


0.5278(77) 


0.02 


0.02 


5- 


16 


1.8(11) 


0.5425(64) 


0.02 


0.025 


5- 


16 


1.7(11) 


0.5567(56) 


0.02 


0.03 


5- 


16 


1.6(11) 


0.5712(50) 


0.02 


0.035 


5- 


16 


1.7(11) 


0.5860(47) 


0.02 


0.04 


5- 


16 


1.9(11) 


0.6010(44) 


0.02 


0.045 


5- 


16 


2.1(11) 


0.6160(41) 


0.02 


0.05 


5- 


16 


2.3(12) 


0.6309(39) 
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TABLE VII: The fitted vector meson mass from the wall-point correlation function for fn S ea, — 

0.03. 



^sca 


val 


fit ran£ 


*e x 2 / dof 


mass 


0.03 


0.01 


6- 


•16 


0.21(32) 


0.537(13) 


0.03 


0.015 


6- 


■16 


0.24(34) 


0.5522(96) 


0.03 


0.02 


6- 


■16 


0.37(42) 


0.5669(79) 


0.03 


0.025 


6- 


■16 


0.55(52) 


0.5809(67) 


0.03 


0.03 


6- 


■16 


0.75(62) 


0.5946(58) 


0.03 


0.035 


6- 


■16 


0.93(70) 


0.6079(52) 


0.03 


0.04 


6- 


■16 


1.07(76) 


0.6213(47) 


0.03 


0.045 


6- 


■16 


1.17(81) 


0.6347(43) 


0.03 


0.05 


6- 


■16 


1.25(85) 


0.6481(40) 



TABLE VIII: The fitted vector meson mass from the wall-point correlation function for ni sea — 0.04. 



m sea , 




fit rai 


ige x 2 /dof 


mass 


0.04 


0.01 


7-16 


1.72(99) 


0.580(25) 


0.04 


0.015 


7-16 


1.8(10) 


0.586(19) 


0.04 


0.02 


7-16 


1.7(10) 


0.591(15) 


0.04 


0.025 


7-16 


1.58(99) 


0.598(11) 


0.04 


0.03 


7-16 


1.51(95) 


0.6084(92) 


0.04 


0.035 


7-16 


1.47(95) 


0.6198(79) 


0.04 


0.04 


7-16 


1.44(97) 


0.6323(70) 


0.04 


0.045 


7-16 


1.41(99) 


0.6455(64) 


0.04 


0.05 


7-16 


1.37(100) 0.6590(59) 
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TABLE IX: The results of a fit to the quark mass dependence of vector meson mass. 

fit X 2 /dof a b c 

dynamical 0.9 (19) 0.448(15) 4.52(47) 

partially quenched 0.34(40) 0.448(13) 1.78(43) 2.78(13) 
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TABLE X: The fitted nucleon mass from the wall-point correlation function. 





fit ran£ 


;e x 2 /dof mass 


0.020 0.010 


6-16 


0.76 (69) 0.679 (19) 


0.020 0.015 


6-16 


0.58 (62) 0.719 (12) 


0.020 0.020 


7-16 


0.90 (80) 0.755 (12) 


0.020 0.025 


7-16 


1.28 (96) 0.789 (10) 


0.020 0.030 


7-16 


1.6 (1.1) 0.8209 (94) 


0.020 0.035 


7-16 


1.8 (1.1) 0.8510 (89) 


0.020 0.040 


7-16 


2.0 (1.1) 0.8796 (85) 


0.020 0.045 


7-16 


2.1 (1.1) 0.9071 (83) 


0.020 0.050 


7-16 


2.2 (1.1) 0.9338 (80) 


0.030 0.010 


7-16 


1.43 (88) 0.747 (30) 


0.030 0.015 


7-16 


1.29 (93) 0.766 (19) 


0.030 0.020 


7-16 


1.13 (84) 0.788 (14) 


0.030 0.025 


8-16 


1.14 (86) 0.816 (16) 


0.030 0.030 


8-16 


1.05 (80) 0.844 (14) 


0.030 0.035 


8-16 


1.01 (77) 0.870 (12) 


0.030 0.040 


8-16 


1.02 (77) 0.896 (11) 


0.030 0.045 


8-16 


1.06 (79) 0.922 (10) 


0.030 0.050 


8-16 


1.12 (81) 0.9481 (95) 


0.040 0.010 


8-15 


1.5 (1.1) 0.748 (44) 


0.040 0.015 


8-16 


0.75 (70) 0.778 (35) 


0.040 0.020 


8-16 


0.53 (56) 0.808 (25) 


0.040 0.025 


8-16 


0.53 (59) 0.835 (19) 


0.040 0.030 


8-16 


0.51 (60) 0.860 (15) 


0.040 0.035 


9-16 


0.33 (54) 0.875 (15) 


0.040 0.040 


8-16 


0.40 (52) 0.910 (11) 


0.040 0.045 


8-16 


0.40 (53) 0.935 (10) 


0.040 0.050 


8-16 


0.47 (59) 0.9602 (92) 



50 



TABLE XI: The ./V* mass from the wall-point correlation function. 





fit range 


X 2 /dof 


mass 


0.020 0.010 


3-5 


0.10 (67) 


1 nnfi 

l.UUD 




0.020 0.015 


5-9 


1.6 (1.4) 


1 D91 




0.020 0.020 


5-9 


1.5 (1.4) 


1 091 


(71 \ 


0.020 0.025 


5-9 


1.4 (1.3) 


1.038 


(59 ) 


0.020 0.030 


5-9 


1.3 (1.3) 


1 061 


(*>()) 


0.020 0.035 


5-9 


1.1 (1.2) 


1 nsfi 

l.Uou 




0.020 0.040 


5-9 


0.9 (1.1) 


1.111 


{61 ) 


0.020 0.045 


5-9 


0.72 (99) 


1.136 


(10\ 
V / 


0.020 0.050 


5-9 


0.54 (87) 


1.161 


(9R\ 

v / 


0.030 0.010 


5-8 


0.6 (1.2) 


1.21 


J 5 ) 


0.030 0.015 


5-8 


0.39 (90) 


0.93 


'1 0\ 


0.030 0.020 


5-9 


1.1 (1.3) 0.950 


{\>o ) 


0.030 0.025 


5-9 


0.8 (1.1) 0.990 


(48 ) 


0.030 0.030 


5-9 


0.64 (98) 1.026 


v ) 


0.030 0.035 


5-9 


0.59 (93) 1.059 


v ) 


0.030 0.040 


5-9 


0.60 (93) 1.089 


v ) 


0.030 0.045 


5-9 


0.62 (94) 1.118 


(IK} 
{ZO ) 


0.030 0.050 


5-9 


0.64 (95) 1.145 


(oi\ 

v J 


0.040 0.010 


5-7 


0.21 (97) 


0.69 


'1 7\ 


0.040 0.015 


5-8 


0.16 (63) 


0.92 


'i n"i 


0.040 0.020 


4-8 


0.19 (52) 


1.020 


I A7\ 


0.040 0.025 


5-9 


0.50 (84) 


1.057 


(49) 


0.040 0.030 


5-9 


0.55 (88) 


1.091 


(39) 


0.040 0.035 


5-9 


0.56 (88) 


1.119 


(33) 


0.040 0.040 


5-9 


0.54 (87) 


1.144 


(28) 


0.040 0.045 


5-9 


0.52 (85) 


1.168 


(25) 


0.040 0.050 


5-9 


0.50 (84) 


1.192 


(23) 



51 



TABLE XII: Physical baryon mass at the light quark mass (fh) by diagonal (m va i = m sea — ► fh) 
or two stage (m va i — > fh, m sea — » fh) extrapolation. The valence extrapolation (m va \ — > fh) results 
are obtained with linear fits using all valence masses in Tables IXl and IXTl 





M N 


X 2 /dof 




X 2 /dof 


fh (diagonal) 


0.605 (26) 


0.5 (1.5) 


0.81 (12) 


1.1 (2.1) 


fh (two stage) 


0.556 (44) 


0.4 (1.3) 


1.00 (13) 


0.7 (1.7) 


0.02 


0.633 (15) 


0.27 (32) 


0.951 (57) 


0.04 (18) 


0.03 


0.687 (21) 


0.024 (67) 0.841 (84) 0.035 (97) 


0.04 


0.704 (35) 


0.05 (11) 


0.892 (73) 


0.50 (48) 



TABLE XIII: Mass ratio of physical iV (N*) and p. and Mn* are calculated with both diagonal 
and two stage extrapolations (Table 1X11)1 . while M p is always from diagonal extrapolation. 

Fit Mn /M p M N */M p 

diagonal 1.329 (59) 1.77 (26) 

two stage 1.221 (98) 2.20 (30) 



52 



TABLE XIV: The fitted pseudo-scalar meson mass from the wall-point correlation function. 



^sca 




fit ran£ 


;e rVdo-f 


mass 


0.02 


0.01 


9- 


•16 


0.79(85) 


0.2160(37) 


0.02 


0.015 


9- 


•16 


1.05(95) 


0.2556(31) 


0.02 


0.02 


9- 


•16 


1.3(11) 


0.2902(28) 


0.02 


0.025 


9- 


•16 


1.6(11) 


0.3213(26) 


0.02 


0.03 


9- 


-16 


1.8(12) 


0.3501(25) 


0.02 


0.035 


9- 


•16 


2.0(13) 


0.3771(24) 


0.02 


0.04 


9- 


•16 


2.2(13) 


0.4026(23) 


0.02 


0.045 


9- 


■16 


2.3(14) 


0.4269(22) 


0.02 


0.05 


9- 


■16 


2.5(14) 


0.4502(21) 


0.03 


0.01 


9- 


•16 


1.6(12) 


0.2240(28) 


0.03 


0.015 


9- 


•16 


1.3(11) 


0.2631(24) 


0.03 


0.02 


9- 


-16 


1.04(94) 


0.2975(22) 


0.03 


0.025 


9- 


•16 


0.90(85) 


0.3287(20) 


0.03 


0.03 


9- 


-16 


0.81(80) 


0.3575(19) 


0.03 


0.035 


9- 


-16 


0.77(77) 


0.3844(18) 


0.03 


0.04 


9- 


-16 


0.76(76) 


0.4098(17) 


0.03 


0.045 


9- 


-16 


0.78(77) 


0.4339(16) 


0.03 


0.05 


9- 


-16 


0.82(79) 


0.4570(15) 


0.04 


0.01 


9- 


-16 


1.06(91) 


0.2254(38) 


0.04 


0.015 


9- 


-16 


1.08(91) 


0.2639(35) 


0.04 


0.02 


9- 


-16 


1.11(91) 


0.2978(33) 


0.04 


0.025 


9- 


-16 


1.14(91) 


0.3287(31) 


0.04 


0.03 


9- 


-16 


1.15(91) 


0.3573(29) 


0.04 


0.035 


9- 


-16 


1.13(91) 


0.3840(27) 


0.04 


0.04 


9- 


•16 


1.09(89) 


0.4094(25) 


0.04 


0.045 


9- 


■16 


1.03(87) 


0.4336(24) 


0.04 


0.05 


9- 


■16 


0.96(84) 


0.4568(23) 
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TABLE XV: The fitted pseudo-scalar meson mass from the axial-vector wall-point correlation 
function. 



w-soa mval n t range x 2 /dof mass 



0.02 


0.01 


9- 


16 


0.62(79) 


0.2154(35) 


0.02 


0.015 


9- 


16 


0.30(56) 


0.2558(28) 


0.02 


0.02 


9- 


16 


0.22(45) 


0.2910(24) 


0.02 


0.025 


9- 


16 


0.27(44) 


0.3227(22) 


0.02 


0.03 


9- 


16 


0.40(50) 


0.3518(21) 


0.02 


0.035 


9- 


16 


0.55(59) 


0.3789(20) 


0.02 


0.04 


9- 


16 


0.72(68) 


0.4045(20) 


0.02 


0.045 


9- 


16 


0.90(76) 


0.4288(19) 


0.02 


0.05 


9- 


16 


1.06(84) 


0.4520(19) 


0.03 


0.01 


9- 


16 


0.56(94) 


0.2221(38) 


0.03 


0.015 


9- 


16 


0.50(77) 


0.2620(32) 


0.03 


0.02 


9- 


16 


0.48(70) 


0.2967(29) 


0.03 


0.025 


9- 


16 


0.48(66) 


0.3280(27) 


0.03 


0.03 


9- 


16 


0.51(65) 


0.3568(25) 


0.03 


0.035 


9- 


16 


0.55(66) 


0.3837(24) 


0.03 


0.04 


9- 


16 


0.61(70) 


0.4091(23) 


0.03 


0.045 


9- 


16 


0.67(74) 


0.4333(22) 


0.03 


0.05 


9- 


16 


0.73(79) 


0.4564(21) 


0.04 


0.01 


9- 


16 


0.40(55) 


0.2241(35) 


0.04 


0.015 


9- 


16 


0.56(64) 


0.2621(29) 


0.04 


0.02 


9- 


16 


0.64(70) 


0.2962(26) 


0.04 


0.025 


9- 


16 


0.62(70) 


0.3273(25) 


0.04 


0.03 


9- 


16 


0.55(68) 


0.3562(23) 


0.04 


0.035 


9- 


16 


0.47(64) 


0.3831(22) 


0.04 


0.04 


9- 


16 


0.40(61) 


0.4086(21) 


0.04 


0.045 


9- 


16 


0.34(60) 


0.4329(20) 


0.04 


0.05 


9- 


16 


0.31(61) 


0.4561(20) 
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TABLE XVI: The pseudo-scalar meson mass and decay constant computed from the pseudo-scalar 
point-point correlation function. 





"1-val 


fit rang 


;e xV dof 


mass 


decay constant 


0.02 


0.01 


9-16 


0.95(81) 


0.2211(28) 


8.80(11) x 10~ 2 


0.02 


0.02 


9-16 


0.39(52) 


0.2938(18) 


9.494(62) x 10~ : 


0.02 


0.03 


9-16 


0.57(63) 


0.3528(19) 


0.10161(78) 


0.02 


0.04 


9-16 


0.54(62) 


0.4051(17) 


0.10720(77) 


0.02 


0.05 


9-16 


0.52(62) 


0.4525(16) 


0.11244(77) 


0.03 


0.02 


9-16 


0.53(66) 


0.3050(26) 


9.740(86) x 10~ : 


0.03 


0.03 


9-16 


0.68(68) 


0.3610(18) 


0.10253(56) 


0.03 


0.04 


9-16 


0.61(69) 


0.4123(20) 


0.10926(76) 


0.04 


0.04 


9-16 


1.4(11) 


0.4087(16) 


0.11059(57) 



TABLE XVII: The pseudo-scalar meson mass and decay constant computed from the axial-vector 
point-point correlation function. 





m va i 


fit rang 


;e x 2 / d °f 


mass 


decay constant 


0.02 


0.01 


9-16 


0.43(54) 


0.2110(76) 


8.97(32) x 10~ 2 


0.02 


0.02 


9-16 


0.47(63) 


0.2891(39) 


9.55(18) x lO- 2 


0.02 


0.03 


9-16 


0.58(76) 


0.3487(40) 


0.1011(20) 


0.02 


0.04 


9-16 


0.73(88) 


0.4014(34) 


0.1058(18) 


0.02 


0.05 


9-16 


0.78(92) 


0.4492(30) 


0.1100(17) 


0.03 


0.02 


9-16 


0.65(69) 


0.3065(59) 


9.75(29) x 10~ 2 


0.03 


0.03 


9-16 


1.5(10) 


0.3610(31) 


0.1032(14) 


0.03 


0.04 


9-16 


0.59(71) 


0.4134(40) 


0.1086(24) 


0.04 


0.04 


9-16 


0.52(64) 


0.4124(25) 


0.1100(16) 
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TABLE XVIII: Results from a fit of the pseudo-scalar meson mass squared to a linear form. 
Only the dynamical, m sea = m va i, points are included. In contrast to the case in the quenched 
approximation, when extrapolated to zero quark mass the result is consistent with zero. This may 
be explained by the absence of a contribution from the quenched chiral logarithm. 



fit ranf 


s;e x 2 /dof c 


b 


Pseudo-scalar 


9-16 


1.0(20) -0.0047(39) x 10~ 3 


4.19(13) 


Axial 


9-16 


0.5(14) -3.1(34) x 10~ 3 


4.12(11) 



TABLE XIX: Parameters from chiral perturbation theory fits to the values of m 2 computed from 
the pseudo-scalar wall-point, and axial- vector wall point. (Tables IXIVI and IXVI respectively) . x 2 
is from uncorrelated fits in mf(= rn seSkiVa \). 



rrif range 


fit ran£ 


re xVdof 2B L 5 -2L 8 


_Zj4 — 2-Lg 


Pseudo-scalar 


0.01-0.03 
0.01-0.04 


9-16 
9-16 


0.12(13) 3.94(27) -1.51(74) x 10~ 4 
1.7(10) 4.18(16) -1.4(44) x 10~ 5 


-1.9(12) x 10~ 4 
-1.17(43) x 10~ 4 


Axial 


0.01-0.03 
0.01-0.04 


9-16 
9-16 


0.22(17) 4.04(28) -1.87(90) x 10~ 4 
1.66(80) 4.23(14) -4.0(48) x 10" 5 


-1.2(11) x 10~ 4 
-7.7(33) x 10~ 5 



TABLE XX: Parameters from next-to-leading order chiral perturbation theory fits to values of f ps 
listed in Table IXVll L, refer to Gasser-Leutwyler low energy constants evaluated at fj, = 1 GeV. 
X 2 is from uncorrelated in mf(= m seaii , m va ]) fits. 



nrif range fit ranj 


;e x 2 / dof a f 


L 5 


u 


0.01-0.03 9-16 
0.01-0.04 9-16 


0.14(32) 5.36(48) x 10" 2 
3.2(18) 4.54(33) x 10~ 2 


7.92(96) x 10~ 4 
7.14(80) x 10~ 4 


7.2(63) x 10~ 5 
1.29(23) x 10~ 4 
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TABLE XXI: Parameters from next-to-next-to-leading order chiral perturbation theory fits to 
values of f ps listed in Table IXVll All quark mass points were used in the fit. Li refer to Gasser- 
Leutwyler low energy constants evaluated at fx = 1 GeV. Ci are 0(p e ) counter-terms, x 2 is from 
uncor related in rrif(= m se a,, rn va \) fits. 

X 2 /dof af L5 U Ci C 2 C12 

0.56(58) 6.43(91) x 10~ 2 5.6(27) x 1(T 4 -2.2(44) x 10~ 4 3.2(76) 1.62(100) 6.8(28) 



TABLE XXII: Parameters from linear fits to values of f ps listed in Table IXVll x 2 is from un- 
corrected in rrif{= m SCSut , m va i) fits. For comparison, results from a two parameter linear fit, 
fp = a f + ci(m + m res ), to the m/ = m sca = m va \ data points are included in the last line. 

nif range fit range x 2 /dof af c\ C2 

0.01-0.04 9-16 0.41(43) 7.81(16) x 10" 2 0.622(20) 0.164(51) 
0.02-0.04 9-16 0.12(72) 7.81(14) x 10" 2 0.783(42) 



TABLE XXIII: The results of fit Eq. [55] without (I = 0) and with (/ ^ 0) the lattice Coulomb 
correction, l6V(r), to data of the static quark potential extracted at t = 5,r £ \r mini r max] = 
[>/3, 8] for 

"wisea — 0.02, 0.03, 0.04. The first and the second errors are the statistical error and the 
estimation of the systematic error due to the selection of the fit ranges, t, [r m i n ,r max ] are swept in 
the region of t = 5, 6, y/2 < r m i n < \/6, and 7 < r max < 9. 





^0 


a a Vq I 


0.02 
0.02 


4.177(22)(99) 0.398(7)(47) 0.0718(11)(49) 0.753(6)(25) (fixed) 
4.126(23) (106) 0.518(18)(57) 0.0665(14)(65) 0.805(10)(38) 0.35(4)(15) 


0.03 
0.03 


4.066(25)(32) 
4.026(26)(62) 


0.368(8)(29) 0.0776(14)(20) 0.728(7)(16) (fixed) 
0.457(21)(55) 0.0736(17)(29) 0.766(11)(26) 0.26(5)(17) 


0.04 
0.04 


4.076(27)(29) 
4.020(29)(34) 


0.399(9)(57) 0.0753(15)(33) 0.749(7)(29) (fixed) 
0.520(23)(28) 0.0699(18)(23) 0.801(12)(13) 0.35(5)(14) 
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TABLE XXIV: The pseudo-scalar B parameter using degenerate valence quarks. 



m sea 


m va i 


t op range 


R (lat) 
-Dps 


0.02 


0.01 


14- 


■17 


0.488(14) 


0.02 


0.02 


14- 


■17 


0.5524(92) 


0.02 


0.03 


14- 


•17 


0.5923(72) 


0.02 


0.04 


14- 


•17 


0.6229(61) 


0.02 


0.05 


14- 


■17 


0.6485(55) 


0.03 


0.01 


14- 


•17 


0.525(14) 


0.03 


0.02 


14- 


•17 


0.5771(83) 


0.03 


0.03 


14- 


■17 


0.6104(64) 


0.03 


0.04 


14- 


•17 


0.6368(55) 


0.03 


0.05 


14- 


•17 


0.6596(50) 


0.04 


0.01 


14- 


•17 


0.512(12) 


0.04 


0.02 


14- 


■17 


0.5747(85) 


0.04 


0.03 


14- 


•17 


0.6133(69) 


0.04 


0.04 


14- 


•17 


0.6425(58) 


0.04 


0.05 


14- 


•17 


0.6662(51) 
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TABLE XXV: The bare pseudo-scalar B parameter using non-degenerate valence quarks, mi and 
771-2, for a dynamical mass of 0.02. 



m sea mi m 2 t op range Bpf^ 

0.02 0.02 0.01 14-17 0.526(12) 

0.02 0.03 0.01 14-17 0.555(11) 

0.02 0.03 0.02 14-17 0.5742(83) 

0.02 0.04 0.01 14-17 0.577(10) 

0.02 0.04 0.02 14-17 0.5929(77) 

0.02 0.04 0.03 14-17 0.6085(67) 

0.02 0.05 0.01 14-17 0.596(11) 

0.02 0.05 0.02 14-17 0.6093(74) 

0.02 0.05 0.03 14-17 0.6231(64) 

0.02 0.05 0.04 14-17 0.6362(58) 



TABLE XXVI: The bare pseudo-scalar B parameter using non-degenerate valence quarks, mi and 
77i2, for a dynamical mass of 0.03. 



m sea mi m 2 t op range Bpf^ 

0.03 0.02 0.01 14-17 0.556(11) 

0.03 0.03 0.01 14-17 0.5799(95) 

0.03 0.03 0.02 14-17 0.5954(73) 

0.03 0.04 0.01 14-17 0.5999(91) 

0.03 0.04 0.02 14-17 0.6116(67) 

0.03 0.04 0.03 14-17 0.6244(59) 

0.03 0.05 0.01 14-17 0.6171(90) 

0.03 0.05 0.02 14-17 0.6262(64) 

0.03 0.05 0.03 14-17 0.6375(56) 

0.03 0.05 0.04 14-17 0.6487(52) 
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TABLE XXVII: The bare pseudo-scalar B parameter using non-degenerate valence quarks, mi and 
m2, for a dynamical mass of 0.04. 



m sea mi m 2 t op range Bpf^ 

0.04 0.02 0.01 14-17 0.551(10) 

0.04 0.03 0.01 14-17 0.5806(97) 

0.04 0.03 0.02 14-17 0.5963(77) 

0.04 0.04 0.01 14-17 0.6045(96) 

0.04 0.04 0.02 14-17 0.6152(74) 

0.04 0.04 0.03 14-17 0.6291(64) 

0.04 0.05 0.01 14-17 0.6248(97) 

0.04 0.05 0.02 14-17 0.6321(72) 

0.04 0.05 0.03 14-17 0.6435(61) 

0.04 0.05 0.04 14-17 0.6551(55) 



TABLE XXVIII: The kaon B parameter from NLO fit, including extrapolation of the sea quark 
mass to the physical point, m sea = rnu g ht- Bo is the value of the B parameter in the chiral limit. 

m sca range m va \ range t op range x 2 /Aoi b b\ - 6 3 b 2 B^ 

0.02-0.03 0.01-0.03 14-17 0.39(39) 0.260(21) 0.527(99) 0.54(26) 0.521(30) 

0.02-0.04 0.01-0.04 14-17 1.27(89) 0.265(11) 0.744(44) 0.26(12) 0.550(16) 

0.02-0.04 0.01-0.05 14-17 2.34(84) 0.2522(86) 0.927(29) 0.27(10) 0.545(14) 

0.02-0.03 0.02-0.03 14-17 0.18(25) 0.258(18) 0.616(64) 0.49(25) 0.525(28) 

0.02-0.04 0.02-0.04 14-17 0.80(93) 0.2594(96) 0.826(33) 0.24(12) 0.547(15) 

0.02-0.04 0.02-0.05 14-17 1.45(72) 0.2476(79) 0.987(25) 0.24(10) 0.543(14) 
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TABLE XXIX: Same as Tabic IXXVHII but including non-degenerate valence quarks. 



m sea range m val range t op range x 2 /dof b b 1 b 2 b 3 B^ 

0.02-0.03 0.01-0.03 14-17 0.20(20) 0.260(22) 0.51(15) 0.56(28) -1(96) x lO" 3 0.524(30) 

0.02-0.04 0.01-0.04 14-17 0.84(74) 0.266(12) 0.713(70) 0.27(13) -1(59) x 10~ 3 0.554(18) 

0.02-0.04 0.01-0.05 14-17 1.36(62) 0.2536(94) 0.853(49) 0.28(11) -3.9(47) x 10~ 2 0.546(16) 

0.02-0.03 0.02-0.03 14-17 9(12) x 10~ 2 0.258(18) 0.25(13) 0.49(25) -0.370(86) 0.498(24) 

0.02-0.04 0.02-0.04 14-17 0.55(79) 0.2591(98) 0.601(54) 0.25(12) -0.220(40) 0.533(15) 

0.02-0.04 0.02-0.05 14-17 0.77(59) 0.2475(83) 0.784(36) 0.25(11) -0.193(31) 0.530(14) 



TABLE XXX: The effective for the plaquette action for an ensemble of 404 quenched DBW2 
configurations versus the observable used in Equation 1731 In each case the operator chosen is 
planar. 

Observable Effective 0p 

1 x 1 9.10039(34) 

2 x 1 7.87358(35) 

2 x 2 6.92080(50) 

3 x 3 6.6815(13) 



TABLE XXXI: The effective, short-distance, gauge action, solved for using the schwinger dyson 
equation with an anzatz for the form of: plaquette (0i), rectangle (0i K 2) and the two hypercubic, 
five-link, loops (05 t i and 05^)- Up to small corrections the results are equal to the bare input 
parameters. 



dynamical mass 


Pixi 


01x2 05,1 


05,2 




0.02 


9.7384(36) 


-1.11392(53) 2.99(56) x 10~ 3 


6.35(77) x 


10~ 3 


0.03 


9.7388(43) 


-1.11523(50) 3.35(67) x 10~ 3 


6.18(93) x 


10~ 3 


0.04 


9.7471(44) 


-1.11531(54) 2.15(73) x 10~ 3 


7.5(10) x 


10~ 3 
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TABLE XXXII: Parameters and results for the stronger coupling dynamical studies 



(3 rrif m p 


X 2 /Dof 


77l res 


# conf. 


0.70 0.026 0.830(5) 


1.4 


0.0094(1) 


32 


0.75 0.022 0.667(7) 


0.8 


0.00405(7) 


41 



TABLE XXXIII: Residual mass versus valence L s for the m se a = 0.04 ensemble 



m val L s = 8 L s = 12 L S = 16 L s = 24 L s = 32 



0.04 4.787(37) x 10" 3 1.347(20) x 10~ 3 5.75(21) x 10" 4 2.13(17) x 10" 4 1.30(15) x 10" 4 
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FIG. 1: Conjugate gradient iteration count for the chronological inverter using the previous 7 
vectors compared with a linear extrapolation of the previous two vectors for the m sea = 0.02 
ensemble. 
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FIG. 2: d(Xff' ,U£') and d max (U}i n) , ') in a trajectory of m sea = 0.02 on 16 3 x 32. The total 
step numbers are 10, 20, 50, 100, 200, and 500. The solid curves are d{ufi\ ujP) while the dashed 
are d max (U^ n ' ,l/jp). 
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FIG. 3: The breaking of reversible dynamics measured by d(U^ ,ujf^) and dmaxiU^^U^) as 
a function of the CG convergence criteria, R CO nv The dotted lines are observed upper bound of 
deviations due to the reunitarization process. The error-bar at R CO nv = le-8 was obtained using 
five configurations. 
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FIG. 4: The residues of CG, Eq. as a function of the number of CG iteration are plotted for 
various numbers of previous solutions, N p , on a typical configuration of the m se a = 0.02 ensemble. 
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FIG. 5: The individual contributions to the total change in the Hamiltonian from the various 
components of the Hamiltonian for the large step-size, old force term simulation described in 
Table Q] The shaded bars represent the trajectories which failed the accept /reject step, while the 
empty bars tally all trajectories. 



67 



3 

a 



c 

o 

1 

ft, 



I 1 I 



I 1 I 



I 1 I 



I 1 I 



_L 



_L 



-1000 -800 -600 -400 -200 200 400 600 800 1000 
1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 



_L 



*.i.;.=.|.=.;f.t=l;.=.;!=.;.i.;i' 



-1000 -800 -600 -400 -200 200 400 600 800 1000 



I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 



_L 



_L 



_L 



_L 



El I "iT. 



_L 



_L 



_L 



-1000 -800 -600 -400 -200 200 400 600 800 1000 



FIG. 6: The individual contributions to the total change in the Hamiltonian from the various 
components of the Hamiltonian for the large step-size, new force term simulation described in 
Table QJ The shaded bars represent the trajectories which failed the accept /reject step, while the 
empty bars tally all trajectories. Note that the magnitude of the differences is dominated by the 
gauge and momentum parts of the Hamiltonian, and the noticeable trend that trajectories for 
which the change in the gauge part of the Hamiltonian is positive are more likely to be rejected. 
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FIG. 7: The history of qq (m va i = mdyn), up to trajectory 1000, for all evolutions. The average 
from trajectory 3000 onwards is shown as a dashed horizontal line. 
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FIG. 8: Plaquette auto-correlation function and integrated auto-correlation length for the m, 
0.02 ensemble. 
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FIG. 9: Auto-correlation function and integrated auto-correlation length for timeslice 12 of the 
correlation function of the time component of the local axial vector current on the m sea = 0.02 
ensemble. 
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FIG. 10: Topological charge history for all the ensembles. Note the correlations of the scale of 
many hundreds of HMC trajectories. 
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FIG. 11: The residual mass plateau plot for m sea = m va i = 0.02. The error-weighted average 
between timeslices 4 and 16 is shown as a horizontal line. 
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FIG. 12: The residual mass extrapolated to m sea — m va i — 0.0. 
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FIG. 13: The renormalisation factor for the local, non-singlet, axial current. This is defined in the 
chiral limit. Here we show the data at finite mass for the fully dynamical points together with the 
results of a linear extrapolation to the chiral limit. 
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FIG. 14: The effective mass in the vector channel for ui sea — Tn va i. 




FIG. 15: The fitted vector mass in the vector channel for m sea = m va i. The solid line denotes a 
linear fit. 
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FIG. 16: The effective mass in the pseudo-scalar channel for fn sea — n^vai — 

0.02 (circles), 0.03 

(squares), and 0.04 (triangles). Here we use the point-point Kuramashi-source correlation function. 
Similar results hold for m sea / m m j. 
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FIG. 17: The effective mass in the pseudo-scalar channel for m sea — n^vai — 

0.02 (circles), 0.03 

(squares), and 0.04 (triangles). Here we use a Coulomb-gauge-fixed wall sourcewith a point sink. 
Similar results hold for m sea ^ m va i. 



79 



0.55 



0.5 



0.45 



t/3 

B 0.4 

> 
• i—i 

£0.35 
pq 



0.3 



0.25 



i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i 



O m =0.02 

w sea 

n m =0.03 

1 sea 

O m =0.04 

sea 



m 



n 2 I— I— I— Lag 

'0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 

timeslice 



FIG. 18: The effective mass in the axial-vector channel for Tn sea — nival — 

0.02 (circles), 0.03 

(squares), and 0.04 (diamonds). Here we use a Coulomb-gauge-fixed wall sourcewith a point sink. 
Similar results hold for m sea ^ m va i- 
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FIG. 19: The pseudo-scalar meson mass, squared, extracted from the wall-point axial-vector cor- 
relation function. m sea = 0.02 (circles), 0.03 (squares), and 0.04 (diamonds). Dashed lines denote 
a next-to-leading order in chiral perturbation theory fit. The dotted line marks the mass of the 
kaon. The solid line shows the results of an extrapolation of the dynamical, and one of the valence, 
quark masses to m. In this case the x-axis represents the remaining valence quark mass. 
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FIG. 20: The pseudo-scalar meson mass squared, extracted from the wall-point axial- vector cor- 
relation function for the dynamical points (m/ = m sca , = m va \), together with the results of a LO 
chiral perturbation theory fit to the quark mass dependence. The experimental value for the square 
of the kaon mass is shown as a dashed line. 
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FIG. 21: The meson decay constant. m sea = 0.02 (circles), 0.03 (squares), and 0.04 (diamonds). 
Dashed lines denote a next-to- leading order in chiral perturbation theory fit, with the solid line 
being the extrapolation to m sea = m. 
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FIG. 22: The meson decay constant. m sea = 0.02 (circles), 0.03 (squares), and 0.04 (diamonds), 
together with the results of a linear fit to the dynamical points. The extrapolated value at m sea = 
—m res is also shown. 
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FIG. 23: The masses of nucleoli and its parity partner for the dynamical points (rn-t — w-sea — ^vai)- 
The solid line is a simple linear fit. Square symbols show the extrapolated results at mj = m. 
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FIG. 24: APE plot of nucleon and its parity partner. Filled symbols show the dynamical DWF 



12| with a ~ 1.3 GeV (squares) 



results, while open symbols show the quenched DWF results 
and a -1 ~ 2 GeV (diamonds). The upper and lower triangles are dynamical DWF results at the 



physical light quark mass (m) with diagonal and two stage chiral extrapolations from Table IXiTll 
Stars indicate experimental values, iV(939) and iV*(1535). 
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FIG. 25: The integrated autocorrelation time, Tint, of the Wilson loops, {W(r, t)), for m sea = 0.02. 
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FIG. 26: The static quark potential of m sea = 0.02 at r, extracted using Eq. |^TJ at t shown as the 
horiozntal axis in the graph. r = v2, 4, 5, and 8. 
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FIG. 27: 5Vl[t) and its corresponding reconstruction from fit data, [y( data )(r) — V CO nt(r)]/l- Data 
extracted on m sea = 0.02 configuration at t = 5, r G [vo, 8] is used in the fit. 
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FIG. 28: The static quark potential extracted from t = 4, 5 and 6 for m sea = 0.02. The black curve 
is the fit result to Eq. EHlwith I = using t = 5 and V3 < r < 8. 
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FIG. 29: £ dependence of the Sommer scale ro for m sea = 0.02,0.03,0.04 as well as the lineary 
extrapolated value at chiral point, m sea = —m res - The fit formula with the lattice Colomb term 
correction .Eg. 1551 with I ^ 0, are used. 
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FIG. 30: Sommer scale ro as a function of m sea and their chiral extrapolation using linear function. 
Error bars are statistical only. 
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FIG. 31: The pseudo-scalar B parameter for m sea = 0.02 and a range of degenerate valence quark 
masses versus the time-slice of the operator insertion. 
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FIG. 32: The pseudo-scalar B parameter for m sea = 0.03 and a range of degenerate valence quark 
masses versus the time-slice of the operator insertion. 
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FIG. 33: The pseudo-scalar B parameter for m sea = 0.04 and a range of degenerate valence quark 
masses versus the time-slice of the operator insertion. 
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FIG. 34: The pseudo-scalar B parameter, extracted from time-slices 14 to 17; m sea ^ m va \. m sea = 
0.02 (circles), 0.03 (squares), and 0.04 (diamonds). The dashed lines are from a fit to Eq. 1651 
evaluated at m sea = 0.04, 0.03, 0.02, and the solid line is an extrapolation, in the dynamical mass, 
to m. The range of quark masses used in this particular fit is 0.02 < mj < 0.04. 
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FIG. 37: Two typical spectral flows: one from the m sea = 0.03 evolution and one from the 
m sea = 0.04 evolution. 
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FIG. 38: A comparison of three spectral flows: all generated with the DBW2 action, using (3 = 0.8, 
(3 = 0.75 and (3 = 0.7 respectively. 



